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ABSTRACT 


The  Affordable  Guided  Airdrop  System  (AGAS)  integrates  a  low-cost  guidance 
and  control  system  into  fielded  cargo  air  delivery  systems.  This  study  evaluated  the 
feasibility  of  this  concept  and  included  the  design  and  execution  of  a  flight  test  program 
to  assess  prototype  system  performance,  as  well  as  modeling  efforts  to  develop  initial 
guidance  and  control  techniques  leading  to  an  evaluation  of  the  feasibility  of  the  AGAS 
concept.  The  flight  test  program  provided  adequate  flight  dynamic  data  for  the  AGAS 
system.  The  wind  measurement  techniques  employed  for  this  effort,  through  the  use  of  a 
"calibration"  parachute  system,  provided  wind  estimates  that  were  not  previously 
available.  Flight  test  data  demonstrated  the  actuator  system  could  provide  glide  ratios  of 
0.4  to  0.5  for  a  flat  circular  parachute.  A  simulation  was  developed  using  a  point  mass 
model  for  parachute  dynamics,  sensor  models,  and  a  Bang-Bang  type  control  system.  Six 
hundred  simulations  demonstrated  that  the  Affordable  Guided  Airdrop  System  shows 
strong  potential  of  providing  a  low-cost  alternative  for  precision  airdrop.  Further  work  is 
recommended  to  implement  a  six-degree  of  freedom  dynamic  model,  assess  the  dynamic 
response  of  the  production  parachute  system,  and  optimize  the  control  algorithms  to 
minimize  fuel  usage. 

This  document  was  prepared  as  a  multimedia  presentation.  There  are  several 
references  to  video  clips  and  animations  throughout  the  printed  document.  Multimedia 
versions  of  this  document  can  be  requested  from  U.S.  Army  Yuma  Proving  Ground, 


STEYP-MT-EA,  Yuma,  A Z  85365. 
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INTRODUCTION 


The  United  States  Air  Force  Science  Advisory  board  was  tasked  to  develop  a 
forecast  of  the  requirements  for  the  most  advanced  air  and  space  ideas  to  project  the  Air 
Force  into  the  next  century.  The  study,  encompassing  all  aspects  of  Air  Force  operations, 
assessed  a  variety  of  technology  developments  critical  to  the  Air  Force  mission.  This 
study  culminated  in  a  report  titled  "New  World  Vistas,  Air  and  Space  Power  for  the  21st 
Century."1  The  study  identified  a  critical  need  to  improve  the  Point-of-Use  Delivery;  that 
is,  getting  the  materiel  where  it  needs  to  be,  when  it  needs  to  be  there.  Airdrop  is  an 
important  aspect  of  Point-of-Use  Delivery.  The  report  indicated  that  immediate 
improvements  are  needed  with  emphasis  provided  by  the  statement:  "In  the  future,  the 
problem  of  airdrop  should  be  treated  as  seriously  as  the  problem  of  bomb  drop." 

To  date,  significant  emphasis  has  been  placed  on  the  development  of  large-scale 
parafoil  systems.  These  systems  provide  the  accuracy  required  with  delivery  from  high 
altitude  and  large  offset  distances.  The  drawback  is  prohibitive  cost  for  each  pound  of 
payload  delivered.  Alternate  approaches  were  required  to  reduce  system  cost.  The  team 
of  the  United  States  Army  and  Air  Force,  The  Boeing  Company,  and  Vertigo, 
Incorporated  is  evaluating  alternative  airdrop  technologies.  These  efforts  include  the 
design  and  development  of  the  Affordable  Guided  Airdrop  System,  which  incorporates  a 
low-cost  guidance,  navigation,  and  control  system  into  fielded  cargo  air  delivery  systems. 
This  study  focused  on  evaluating  the  feasibility  of  the  AGAS  concept  and  encompassed 
the  design  and  execution  of  a  flight  test  program  to  assess  dynamic  response  of  a  flat 
circular  parachute,  to  the  design  of  initial  guidance  and  control  techniques,  and  to 
evaluate  the  feasibility  of  the  AGAS  concept. 
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H.  SYSTEM  DESCRIPTION 

A.  AFFORDABLE  GUIDED  AIRDROP  SYSTEM  (AGAS) 

The  Affordable  Guided  Airdrop  System2  (AGAS)  is  being  evaluated  as  a  low-cost 
alternative  for  meeting  the  military's  requirements  for  precision  airdrop.  Designed  to 
bridge  the  gap  between  expensive  high  glide  ratio  parafoil  systems  and  uncontrolled 
(ballistic)  round  parachutes,  the  AGAS  concept  offers  the  benefits  of  high  altitude 
parachute  releases  but  cannot  provide  the  same  level  of  offset  from  the  desired  impact 
point  (IP)  as  high-glide  systems.  The  design  goal  of  the  AGAS  development  is  to 
provide  a  Guidance,  Navigation,  and  Control  (GNC)  system  that  can  be  placed  in-line 
with  existing  fielded  cargo  parachute  systems  (G-12  and  G-ll)  and  standard  delivery 
containers  (A-22).  The  system  is  required  to  provide  an  accuracy  of  100  meters,  Circular 
Error  Probable  (CEP),  with  a  desired  goal  of  50  meters  CEP.  No  changes  to  the 
parachute  or  cargo  system  are  allowed. 

The  current  design  concept  includes  implementation  of  a  commercial  Global 
Positioning  System  (GPS)  receiver  and  a  heading  reference  as  the  navigation  sensors,  a 
guidance  computer  to  determine  and  activate  the  desired  control  input,  and  the 
application  of  Pneumatic  Muscle  Actuators  (PMAs)  to  effect  the  control.  The  navigation 
system  and  guidance  computer  would  be  secured  to  the  existing  container  delivery 
system  while  the  PMAs  would  be  attached  to  each  of  four  parachute  risers  and  to  the 
container.  Figure  1  illustrates  the  concept.  Control  is  affected  by  lengthening  a  single  or 
two  adjacent  actuators.  The  parachute  deforms  creating  an  asymmetrical  shape, 
essentially  shifting  the  center  of  pressure,  and  providing  a  drive  or  slip  condition.  Upon 
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deployment  of  the  system  from  the  aircraft,  the  guidance  computer  would  steer  the 
system  to  a  pre-planned  trajectory.  This  concept  relies  on  the  ability  of  sufficient  drive  to 
be  produced  to  overcome  errors  in  wind  estimation  and  the  point  of  release  of  the  system 
from  the  aircraft. 


B.  PARACHUTE 

The  C-9  parachute  was  selected  for  this  feasibility  demonstration  due  to  its 
availability  and  representation  of  the  larger  cargo-type  parachutes  (G-l  1  and  G-12)  on 
which  this  system  will  ultimately  be  used.  Although  the  C-9  was  initially  designed  as  an 
ejection  seat  parachute,  it  is  a  standard  flat  circular  parachute  as  are  the  larger  G-l  1  and 
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G-12  cargo  parachutes.  A  flat  circular  parachute  is  one  that  when  laid  out  on  the  ground 
forms  a  circle. 


The  reference  diameters  of  these  chutes  are  28  feet  (C-9),  64  feet  (G-12),  and  100 
feet  (G-l  1).  The  reference  area  of  the  C-9  parachute  is  taken  to  be  the  total  surface  area 
of  the  canopy  (a  circle  of  28  foot  diameter)  and  is  615.8  square  feet.  The  C-9  is  static¬ 
line  deployed  and  utilizes  28  suspension  lines  connecting  to  four  risers. 

A  cargo  box  was  suspended  from  the  system  and  housed  the  remote  control 
system,  control  actuators,  and  instrumentation  system. 

C.  ACTUATORS 

Vertigo,  Incorporated  developed  Pneumatic  Muscle  Actuators3  (PMAs)  to  effect 
the  control  inputs  for  this  system.  The  PMAs  are  braided  fiber  tubes  with  neoprene  inner 
sleeves  that  can  be  pressurized.  Uninflated  PMAs  as  installed  on  a  scaled  system  are 
shown  in  Figure  3 .  Upon  pressurization,  the  PMAs  contract  in  length  and  expand  in 
diameter. 
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Figure  3.  Pneumatic  Muscle  Actuators  (PMAs) 

With  four  independently  controlled  actuators,  two  of  which  can  be  activated 
simultaneously,  eight  different  control  inputs  can  be  affected.  For  this  demonstration,  a 
throw  of  approximately  3  feet  was  selected  [measured  to  be  3.25  feet  during  ground 
testing].  When  depressurized,  the  PMAs  are  completely  flexible  allowing  for  efficient 
packing  of  the  actuators  with  the  parachute. 


Figure  4.  Packing  the  Parachute  and  Actuators 


The  concept  employed  for  the  AGAS  is  to  fully  pressurize  all  actuators  upon 
successful  deployment  of  the  parachute.  To  affect  control  of  the  system,  one  or  two 
actuators  are  depressurized.  This  action  "deforms"  the  parachute  creating  drive  in  the 
opposite  direction  of  the  control  action.  Figure  5  illustrates  this  action  while  the  Figure  6 
illustrates  the  parachute  deformation  upon  control  actuation. 


Figure  6.  Parachute  Response  to  Control  Input  (video  clip) 
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D.  MASS  AND  CENTER  OF  GRAVITY 

The  weights  of  the  major  system  components  were  determined  through  direct 
measurement  or  calculation  using  known  material  weights  (e.g.  parachute  canopy  fabric 
weight  of  1.1  ounces  per  square  yard).  The  major  system  component  weights  are 
summarized  in  Table  1. 

The  theoretical  center  of  gravity  of  the  system  was  found  to  be  1 1.4  feet  vertically 
upward  from  the  center  of  the  payload  for  a  standard  sea-level  atmosphere.  This 
calculation  includes  the  mass  of  the  trapped  air  in  the  canopy  and  as  such  is  dependent  on 
the  air  density.  The  relationship  of  center  of  gravity  to  air  density  is  as  follows: 

(V  4,  =  VP(pVg  +  Wp)  +  laWa  +ltWt  +lcWc]/W^m 

Substituting  known  parameters  equation  1  reduces  to: 


"o' 

leg  =  4323/0  +  0.5  feet  and  Pa  =  0 

l 

L  **  J 

0 

=  0 

4323/7 +  .5 

Component 

Weight  (pounds) 
sea  level  -  standard  atmosphere 

Parachute 

7.5 

Suspension  Lines 

4.5 

Actuators  (total  for  four) 

13.3 

Payload 
cargo  box 
instrumentation 
actuator  valves/nitrogen  tanks 

320.7 

Total  System 

346.0 

Table  1.  System  Weight 


E.  MOMENTS  OF  INERTIA 

The  moments  of  inertia  of  the  parachute  (including  trapped  air  mass),  payload, 
suspension  lines,  and  actuators  were  determined  as  shown  below.  A  summary  of  these 
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values  is  shown  in  Table  2.  Due  to  symmetry  of  the  parachute  and  payload  system,  all 
cross  products  of  inertia  where  assumed  to  be  zero. 


Component 

Moment  of  Inertia  (slugs-ff2) 
sea  level  -  standard  atmosphere 

Ixx 

In 

hi 

Parachute 

1865.9 

1865.9" 

211.3 

Suspension  Lines 

neglected 

Actuators  (total  for  four) 

neglected 

Payload 

1328.3 

1328.3 

39.9 

cargo  box 

instrumentation 

actuator  valves/nitrogen  tanks 

Total  System 

3194.2 

3194.2 

251.2 

Table  2.  Moments  of  Inertia 


1.  Parachute 

When  inflated,  a  flat  circular  parachute  approximates  a  hemisphere  with  a  radius 
(r).  The  moments  of  inertia  can  be  found  as  follows: 

2 

(2)  I XX  ~^YY  ~^2Z  =~^mr  ' 

These  values  then  must  be  translated  to  the  principal  axis  located  at  the  center  of 
gravity  at  the  system.  The  parallel  axis  theorem4  is  applied  as  follows: 

(V  I xx-  =  In.  =Ixx+  ™dx2  and  1^.  =  I 

2  2 

(4)  I xx-  ~  4 yy’  =  ~mr  +  mdx  where:  dx=dp-lcg. 

The  moment  of  inertia,  like  the  center  of  gravity,  varies  with  air  density  as  the 
trapped  air  mass  is  a  significant  contributer  to  this  term.  Substituting  m  =  mp+  pVg 
(where  mp  is  the  mass  of  the  parachute  material  and  all  known  parameters: 

I  xx-  =  Ijj-  =  734436/?  + 1449.0 
I ^  =83152/9  +  53.5 
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2.  Payload 

The  cargo  box  dimensions  were  2  feet  by  2  feet  by  2  feet.  It  is  assumed  that  the 
center  of  gravity  of  this  box  was  located  at  the  center  of  the  box.  The  moments  of  inertia 
can  be  found  as  follows: 

(V  I xx  =  In  =  Izz  +/2);  where  a  =  l . 

Again  applying  the  parallel  axis  theorem,  the  moments  of  inertia  about  the 
principal  axis  were  found.  The  parallel  axis  theorem  is  applied  as  follows: 

(6)  I xx  I n‘  ~  I xx  ntdx  and  1 77.  —  1 ^ . 

3.  Actuators  and  Suspension  Lines 

The  actuators  and  suspension  lines  were  treated  as  slender  rods  with  the  body  axis 
located  at  the  end  of  the  rod.  The  moments  of  inertia  were  found  as  follows: 

(V  I XX.  =  In=  and  4z  =  °- 

These  values  for  this  system  were  found  to  be  extremely  small  as  compared  to  the 
values  for  the  parachute  and  payload  and  therefore  were  considered  negligible. 
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m.  CONTROL  SYSTEM 


For  an  airdrop  mission,  the  aircrew  will  determine  the  Computed  Air  Release 
Point  (CARP)  based  on  the  best  wind  estimate  available  at  that  time.  The  aircraft  will 
then  be  navigated  to  that  point  for  air  delivery  of  the  materiel.  Should  the  wind  estimate 
and  calculation  of  the  predicted  release  point  be  perfect  and  the  aircrew  gets  the  aircraft 
to  the  precise  release  point,  then  the  parachute  would  fly  precisely  to  the  target  without 
control  inputs.  However,  wind  estimation  is  far  from  a  precise  science.  The  calculation 
of  the  CARP  relies  on  less  than  perfect  estimates  of  parachute  aerodynamics  and  the 
flight  crews  cannot  possibly  precisely  hit  the  predicted  release  point  for  each  airdrop 
mission.  Therefore,  the  AGAS  control  system  design  must  help  overcome  these  potential 
errors. 

A  glide  ratio  (ratio  of  horizontal  to  vertical  velocity)  demonstrated  in  flight  test 
was  approximately  0.4-0.5.  Considering  this  relatively  low  glide  ratio  and  a  descent  rate 
of  approximately  25  feet  per  second,  the  AGAS  can  overcome  only  a  twelve  foot  per 
second  (approximately  7  knots)  wind.  It  is  therefore  imperative  to  implement  the  system 
to  overcome  poor  estimates  in  the  wind  and  not  try  to  steer  the  system  against  the  entire 
wind.  In  other  words,  the  drive  of  the  system  is  insufficient  to  attempt  to  fly  straight  to 
the  target  but  is  sufficient  to  overcome  up  to  a  twelve  foot  per  second  error  in  the  wind 
estimate.  For  this  reason,  a  trajectory  control  approach  was  selected. 

A  pre-planned  trajectory,  based  on  the  best  wind  estimate  available,  must  be 
determined  and  provided  to  the  guidance  computer.  The  GPS  navigation  system  will 
provide  a  continuous  position  of  the  system.  The  guidance  computer  will  compare  the 
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actual  horizontal  position,  at  the  system's  current  altitude,  to  the  planned  trajectory  at  that 
altitude.  This  represents  the  position  error  (Pe)  at  the  current  time. 

A  tolerance  cone  is  established  about  the  planned  trajectory  (Figure  7)  starting  at 
600  feet  at  the  beginning  of  the  trajectory  and  gradually  decreasing  to  60  feet  at  ground 
level.  Should  the  position  error  be  outside  this  tolerance,  a  control  is  activated  to  slip  the 
system  back  to  the  planned  trajectory.  When  the  system  is  within  30  feet  of  the  planned 
trajectory  the  control  is  disabled  and  the  parachute  drifts  with  the  wind.  Thirty  feet  was 
selected  to  encompass  approximately  1-sigma  of  the  GPS  errors  (Selective  Availability 
off). 


Figure  7.  Control  Concept  (double  click  to  see  animation) 


As  outlined  above,  the  control  system  relies  on  the  current  horizontal  position 
error  to  determine  if  control  input  is  required.  This  position  error  (Pe)  is  determined  in 
inertial  space  and  is  then  rotated  to  the  body  axis  using  an  Euler  angle  rotation  with 
heading  only.  The  resultant  body-axis  error  (Pb)  is  then  used  to  identify  which  control 
input  must  be  activated  as  shown  in  Equation  8.  Two  components  are  returned,  a  +  or  - 
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for  the  x-axis  and  a  +  or  -  for  the  y-axis.  It  was  assumed  for  this  simulation  that  a  +x 
would  activate  control  A,  a  -x  activates  control  C,  a  +y  activates  control  B,  while  a  -y 
activates  control  D.  The  actual  rigging  of  the  operational  system  must  align  these  control 
actuators  to  the  compass  reference  line  to  ensure  proper  control.  Control  A  is  assumed  to 
be  aligned  with  the  compass  zero  reference  line. 

(8)  input  =  sign 

The  magnitude  of  this  calculation  is  used  to  determine  if  the  selected  control  input 
will  be  activated.  If  the  magnitude  is  greater  than  0.3,  then  that  control  is  activated.  This 
will  allow  the  activation  of  a  single  control  input  or  two  simultaneous  control  inputs: 
Preventing  activation  of  control  when  this  value  is  less  than  0.3  will  ensure  that 
unnecessary  control  inputs  are  not  activated  when  the  predominate  error  is  in  a  single 
direction.  Figure  8  illustrates  the  region  of  active  control. 


Figure  8.  Control  Activation 
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IV.  TESTING 


The  objectives  of  the  test  program  were  to  investigate  the  apparent  mass  effects 
on  flat  circular  parachutes  and  to  obtain  sufficient  data  to  assess  the  feasibility  of  the 
AGAS  concept.  To  meet  these  objectives,  ground  testing  of  the  C-9  parachute  and  flight 
testing  of  the  prototype  AGAS,  using  remotely  controlled  activation  of  the  actuators, 
were  conducted. 

A.  INVESTIGATION  OF  APPARENT  MASS  EFFECTS 

To  illustrate  the  effects  of  apparent  mass  on  the  dynamics  of  a  parachute,  a  test 
program  was  established  to  collect  the  forces  along  the  z-axis  of  the  parachute  during 
acceleration.  The  concept  defined  by  Vertigo,  Incorporated5  to  evaluate  the  apparent 
mass  terms  was  implemented.  The  C-9  parachute  was  attached  to  a  tower  installed  on  the 
vehicle  with  two  attaching  risers  and  towed  behind  a  ground  vehicle  (Figure  9). 


The  vehicle  first  accelerated  to  a  speed  of  approximately  10  feet  per  second. 
After  a  brief  period  of  constant  velocity,  the  vehicle  then  accelerated  at  a  near  constant 
acceleration.  Figure  10  shows  a  typical  velocity  profile  for  these  tests. 
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Figure  10.  Tow  Test  Velocity  Profile 


The  forces  in  the  attaching  risers  were  measured  with  strain  gauges.  In  addition,  a 
"box"  was  constructed  on  top  of  the  tower  to  measure  the  vertical  and  lateral  forces  on 
the  riser  attachment  points.  By  measuring  these  forces,  an  estimate  of  the  parachute 
angle  of  incidence  could  be  derived.  Figure  1 1  illustrates  the  instrumentation 
configuration  at  one  riser  attachment  point.  The  velocity  was  measured  with  differential 
GPS  and  an  anemometer. 

To  estimate  the  actual  apparent  mass  parameters,  the  motion  of  the  parachute 
would  be  ideally  fixed  in  one  direction.  In  the  Vertigo  test,  the  parachute  used  was  very 
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stable.  Therefore,  the  assumption  of  motion  in  only  one  axis  was  very  reasonable.  With 
the  C-9  parachute  significant  oscillations  were  observed. 


Figure  11.  Instrumentation  Configuration:  Riser  Attachment 
A  plot  of  the  calculated  incidence  angle  (<j>)  demonstrates  the  oscillatory  nature  of 
this  parachute.  Therefore,  estimation  of  the  apparent  mass  coefficients  using  these 
techniques  was  not  possible. 


Figure  12.  Incidence  Angle  (<j>) 
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Even  though  precise  apparent  mass  coefficients  cannot  be  found,  insight  to  the 
apparent  mass  contributions  can  be  gained  by  analyzing  these  data.  Figure  13 
demonstrates  the  influence  of  apparent  mass  on  parachute  dynamics.  The  total  force 
along  the  parachutes  z-axis,  as  measured  in  the  risers,  is  plotted  along  with  the  steady 
state  drag  calculations  as  found  by: 

(9)  As  =  —  pV2  • S0-CD ;  where  S0  is  the  reference  area  of  the  parachute  and  C&  is 
the  steady-state  drag  coefficient  (0.68). 


Figure  13.  Apparent  Mass  Effects  on  Parachute  Dynamics 


The  figure  shows  increased  measured  force  above  the  predicted  steady-state  drag 
along  the  parachute  when  the  parachute  is  accelerating.  Note  during  the  first 
accelerations  that  the  measured  force  did  not  rise  above  the  steady  state  drag  estimate. 
This  is  likely  due  to  the  parachute  inflation  and  the  greatly  changing  angle  of  incidence 
during  this  inflation.  Concentrating  on  the  area  where  the  parachute  has  settled  into  pure 
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oscillations,  the  apparent  mass  effects  can  be  observed.  Figure  14  presents  an  excerpt  of 
data  for  this  time: 


These  data  demonstrate  that  the  apparent  mass  effects  must  be  considered  in  flight 
dynamic  modeling.  Even  with  the  relatively  small  accelerations  experienced  here,  the 
effect  is  very  large  and  varies  greatly  with  acceleration  of  the  parachute  (as  observed  with 
the  oscillating  parachute).  Techniques  for  estimating  apparent  mass  coefficients  are 
presented  in  the  system  modeling  section  of  this  report. 

B.  FLIGHT  TESTING 

The  flight  test  effort  focused  on  the  collection  of  flight  dynamic  data  to  support 
modeling  of  the  AGAS  concept.  The  flight  test  effort  was  conducted  with  four  actuators 
in-line  with  a  C-9  parachute  and  a  one-half  scale  container  delivery  system.  The 
actuators  were  activated  using  a  manual  radio  control  system.  Flight  dynamic  data  were 
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obtained  including  the  position,  velocity,  acceleration,  attitude,  and  attitude  rates  of  the 
system.  It  was  necessary  to  correlate  these  data  with  control  inputs.  Therefore,  the  state 
of  control  activation  was  monitored.  Parachute  performance  is  significantly  influenced 
by  the  winds.  It  was  critical  to  this  effort  to  measure  the  winds  as  precisely  as  possible. 

C.  INSTRUMENTATION 

Ideally,  both  the  parachute  and  payload  would  have  been  instrumented  to  collect 
all  necessary  data.  However,  the  state-of-the-art  in  instrumentation  is  not  yet  sufficient  to 
adequately  instrument  the  parachute  itself.  As  a  result,  only  the  payload  could  be 
instrumented.  A  custom  instrumentation  system  was  developed  and  it  included  a 
differential  GPS  system  for  precise  position  and  velocity,  3-axis  accelerometers  for 
acceleration,  and  an  Attitude  Heading  Reference  System  (AHRS)  for  3-axis  attitudes  and 
attitude  rates.  Pressure  transducers  were  put  in  line  with  the  pneumatic  actuators  to 
monitor  their  action.  Figure  15  illustrates  the  instrumentation  design.  A  summary  of  the 
major  instrumentation  components  is  provided  in  Appendix  D. 


Figure  15.  Instrumentation  Block  Diagram 


20 


The  critical  part  of  this  design  is  the  time  synchronization  of  data  from  all 
sources.  To  achieve  this,  the  IRIG  time  generator  was  synchronized  with  GPS  time  using 
the  Havequick  Time  Interface  and  the  1-Pulse  Per  Second  time  sync.  The  AHRS  data 
were  available  at  its  interface  approximately  71  times  per  second.  This  is  the  fastest  rate 
for  all  the  sensors.  Therefore,  these  data  were  used  as  the  key  for  capturing  data  from  all 
sensors.  When  an  AHRS  message  was  first  detected  by  CPU1,  the  first  record  was  time- 
tagged  and  the  entire  message  was  sent  to  the  recorder.  At  that  time,  the  Analog-Digital 
(A-D)  converter  was  polled  for  its  data.  The  first  record  from  the  A-D  converter  was 
time-tagged  and  the  message  sent  to  the  data  recorder.  The  GPS  data  from  the 
differential  receiver  were  independently  captured  (with  an  embedded  time-tagging  using 
GPS  time)  and  recorded  on  a  separate  PC-Card.  The  synchronization  of  the  data  was 
validated  in  two  ways.  First,  it  was  desired  to  have  a  discrete  event  that  would  effect  all 
data  sensors.  The  obvious  event  was  ground  impact,  which  resulted  in  immediate 
changes  in  the  data  from  each  sensor.  This  discrete  event  showed  that  the  data  from  the 
independent  sensors  were  time  synchronized  to  less  than  200  milliseconds  (the  rate  of 
GPS  data).  Next,  the  acceleration  data  were  integrated  to  estimate  velocity.  This 
estimated  velocity  was  then  compared  to  the  GPS  velocity  (after  rotation  from  the  body 
reference  frame  to  the  inertial  reference  frame  using  the  AHRS  attitude  data).  These 
results  also  showed  the  data  were  time-tagged  to  within  200  milliseconds. 

D.  WIND  ESTIMATION 

Two  methods  of  wind  estimation  are  presented.  The  first  includes  the  accepted 
standard  of  the  Radiosonde  Wind  Measuring  System  (RAWIN)  system  used  throughout 
the  test  community.  A  RAWIN  balloon  was  launched  at  approximately  one-hour 
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intervals  near  the  release  time  in  the  vicinity  of  the  Drop  Zone.  This  system  provided  a 
direct  measurement  of  the  winds  as  a  function  of  altitude.  Although  an  accepted 
standard,  the  RAWIN  system  has  limitations  in  airdrop  operations.  The  largest  problem 
is  that  real-time  winds  are  not  available.  The  balloon  must  be  launched  and  data 
processed  resulting  in  approximately  a  one-hour  delay.  Figure  16  illustrates  the 
magnitude  of  wind  changes  over  time.  Data  from  three  RAWIN  launches  are  presented 
as  a  function  of  altitude.  It  is  clear  that  as  the  parachute  gets  close  to  the  ground,  the 
wind  changes  can  be  significant. 


Figure  16.  Wind  Changes  Over  Time 


A  second  method  of  deploying  a  "calibration  system"  just  prior  to  release  of  the 
test  payload  was  implemented.  A  tri-lobe  parachute  was  used  with  a  reference  drag  area 
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(CdS)  of  90.88.  Oscillations  of  this  parachute  were  observed  to  be  very  small.  The 
system  was  weighted  at  58  pounds  to  provide  approximately  the  same  descent  rate  as  the 
C-9  system.  The  calibration  system  was  instrumented  with  differential  GPS  on  the 
payload.  Initially,  the  wind  estimate  was  simply  taken  as  the  ground  velocity  as 
measured  by  the  GPS.  This  approach  effectively  considers  the  calibration  system 
massless  and  therefore  does  not  account  for  changes  in  momentum. 

To  validate  the  ability  to  simply  use  the  measured  ground  track  velocity  as  the 
wind  estimate,  a  model  of  the  calibration  parachute  system  was  developed.  A  point-mass 
system  was  assumed  with  the  only  forces  on  the  system  being  drag  and  weight.  The 
applicable- equations  of  motion  are: 


(10)  X  =  (m  +  an)u  = -D cosy  cosy/ 

(11)  7  =  (m  +  )v  =  -D  cosy  sin  y/ 

(12)  Z  =  (m  +  an)w  =  -Dsmy  +  W 

where:  W  is  the  calibration  system  weight,  D  is  drag  (D  =  qCDS ),  y  and  y/  are 
the  flight  path  angle  and  yaw  angle  respectfully.  The  reference  angles  are: 


Substituting,  rearranging  terms,  and  putting  in  state  space  form: 


where:  are  the  apparent  mass  terms,  here  assumed  to  be  constant 
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Assuming  no  rotation  between  the  fixed  earth  reference  and  the  system's  body 
axis,  the  ground  velocity  can  be  determined. 

Vg  ~VAJrVw\  where  Vg  ,  VA,  and  Vw  are  the  ground  velocity,  parachute  velocity 
relative  to  the  local  air  mass  (airspeed),  and  wind  velocity,  respectively. 

(15)  VG=VA+VW 
Applying  15  to  14: 

u  m+an  0  0  T1  f  r«l  Tot 

(16)  VG  -  v  =  0  m  +  a22  0  v  +  0  >  +  Vw 

*io  0  0  m  +  a33J  T  w  W 

Equation  16  can  be  solved  numerically  to  estimate  the  system  response  to  the 
estimated  winds.  This  estimate  can  then  be  compared  to  the  measured  velocity  of  the 
system.  Using  the  measured  ground  track  velocity  as  the  initial  wind  estimate,  the 
modeled  ground  track  was  determined.  The  difference  between  this  modeled  ground 
track  and  the  actual  measured  ground  track  reflects  the  errors  in  the  wind  estimate.  The 
magnitude  of  these  errors  indicates  the  significance  of  accounting  for  momentum  changes 
caused  by  changes  in  the  wind.  A  Simulink®  model  was  implemented  (Figure  17). 


Figure  17.  Simulink®  Realization  for  Wind  Estimation 


®  Registered  trademark  of  Mathworks,  Inc.,  Natick,  MA 
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Figure  18  presents  the  results  of  this  validation.  Use  of  the  measured  ground 
track  velocity  as  the  wind  estimate  resulted  in  errors  of  less  than  0.3  feet  per  second. 


Figure  18.  Wind  Estimation  Results 


To  further  assess  this  technique,  the  equations  for  wind  estimation  utilized  with 
hurricane  dropsondes6  was  applied.  In  this  application,  the  3-dimensional  velocity  of  a 
dropsonde  is  determined.  The  horizontal  velocity  components  of  the  air  mass  are  then 
found  as  follows: 

.  xz  .  yz 

(17)  xw^x - andyw»y - ; 

g  g 

To  apply  equation  17,  a  Simulink®  model  was  implemented  (Figure  19).  The 
correction  that  would  be  applied  to  the  horizontal  velocity  is  presented  in  Figure  20. 
These  data  show  that  the  corrected  wind  estimate,  using  the  hurricane  approach,  differs 
from  the  measured  ground  track  by  less  than  0.3  feet  per  second.  More  errors  are  seen 
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close  to  ground  impact  where  additional  shears  are  present.  This  is  not  of  consequence 
for  this  study  as  the  flight  dynamic  data  of  interest  is  at  much  higher  altitude. 


Figure  19.  Hurricane  Wind  Estimation  Implementation 


Figure  20.  Hurricane  Wind  Estimation  Results 
These  results  indicate  that  the  momentum  effects  can  be  ignored  for  wind 
estimation  for  the  selected  parachute  system.  Figure  21  illustrates  the  comparison  of  the 
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wind  estimate  to  the  winds  measured  by  the  RAWIN.  Recall  that  the  RAWIN  balloon 
was  launched  only  every  hour.  The  closest  RAWIN  data  were  used  for  this  comparison. 


These  results  demonstrate  this  technique  will  provide  significantly  better 
estimates  of  winds  than  using  the  RAWIN  system.  By  adjusting  the  weight  of  the 
calibration  system  to  match  the  descent  rate  of  the  test  item,  the  two  parachutes  will  be 
subjected  to  the  same  (as  close  as  possible)  atmospheric  conditions.  Using  the  measured 
GPS  ground  track  velocities  is  an  adequate  approximation  for  wind  estimation.  Other 
techniques,  such  as  that  presented  above,  may  provide  some  refinement  on  the  wind 
estimates,  but  the  difference  is  likely  to  be  insignificant  for  most  testing.  The  key  to 
application  of  this  technique  is  the  use  of  a  very  stable  parachute  due  to  the  reductions  in 
apparent  mass  effects  resulting  from  oscillations. 
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E.  DATA  REDUCTION 


Data  from  all  sources  were  time  correlated  using  the  AHRS  time  tag.  As  these 
data  were  recorded  at  the  fastest  rate,  all  other  data  parameters  were  correlated  with  this 
time  field.  A  linear  interpolation  was  used  to  estimate  the  acceleration,  actuator 
movement,  and  trajectory  of  the  AGAS  system  at  the  time  of  the  AHRS  measurement. 
The  data  recorded  on  the  calibration  parachute  (differential  GPS  only)  were  correlated  to 
the  system  data  again  using  linear  interpolation.  Instead  of  time  correlating,  the 
calibration  data  were  correlated  to  the  AGAS  trajectory  data  with  the  AGAS  altitude.  As 
previously  discussed,  the  velocity  data  from  the  calibration  parachute  were  taken  as  the 
wind  without  any  processing  (excluding  the  vertical).  For  the  vertical  winds,  the  mean 
and  density  effects  were  removed  from  the  calibration  parachute  vertical  velocity. 

The  performance  of  the  C-9  parachute  is  characterized  by  significant  oscillations 
(Figure  22)  in  both  pitch  and  roll.  Due  to  the  increased  mass  of  the  trapped  air,  the  larger 
G-12  parachute  system  is  not  expected  to  oscillate  as  severely  as  this  parachute. 


Figure  22.  Parachute  Oscillations  (video  clip) 
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The  system  instrumentation  was  all  located  in  the  payload  and  therefore  measured 
the  movement  of  the  payload  induced  from  the  oscillations.  The  velocity  data  best 
demonstrates  this  character  (Figure  23). 


Figure  23.  Data  Characteristics  Due  to  Oscillations 


To  better  represent  the  operational  system  during  this  feasibility  test,  attempts 
were  made  to  modify  the  C-9  parachute  to  minimize  the  oscillations  by  cutting  large 
symmetrical  slots  in  the  canopy.  This  modification  proved  very  stable  during  tow  testing 
as  well  as  the  initial  test  drop.  The  system  was  then  dropped  with  one  riser  extended 
representing  a  single  control  input.  The  objective  of  this  drop  was  to  qualitatively  assess 
the  amount  of  drive  that  could  be  obtained  with  this  modified  parachute.  There  was 
insufficient  indication  that  this  modification  allowed  enough  drive  for  the  guided  system. 
Had  the  results  been  favorable,  additional  instrumented  drops  with  the  actuator/control 
system  would  have  been  accomplished. 


29 


Since  the  trials  with  the  modified  parachute  did  not  present  the  desired  results,  the 
oscillation  in  the  data  needed  to  be  dealt  with  prior  to  parameter  estimation.  These  data 
must  be  corrected  to  the  location  of  the  system's  center-of-gravity. 

Starting  with  the  position  corrections: 

(18)  Pcc  -  Pl +UbPL  ’  where  Pl  is  the  measured  position  of  the  payload  and  L  is  the 

lever-arm  from  the  location  of  the  instrumentation  to  the 
center  of  gravity  in  body-axis  coordinates  {B } . 

The  velocity  correction  is  then  given  by: 

(19)  Vco  =  4,  =  P^RL+^L 

at 

Noting:  BR=BRS(a> ) 

(20)  ubRL  =  vbRS  (to  )L  =  BR(a>xL  ) 

Substituting  (20)  into  (19): 

(21)  VCG  =  PL  +ubR{coxL)+ubR  ~  L 

at 

j 

For  rigid  bodies,  L  is  assumed  constant  and  therefore  —  L  =  0 

dt 

Applying  this  condition  to  a  parachute  did  not  yield  adequate  results,  as  the 
oscillations  were  still  apparent  in  the  velocity  data.  It  is  presumed  that  in  an  oscillating 
parachute,  the  forces  on  the  parachute  are  changing  due  to  the  linear  acceleration  induced 
by  the  oscillations  (apparent  mass  effects)  while  the  forces  on  the  payload  are  not 
varying.  With  this  imbalance  of  forces  in  the  parachute  and  payload,  an  effective  change 
in  center-of-gravity  of  the  system  results. 

To  illustrate  this  situation,  the  lever-arm  correction  is  determined  by  rearranging 
Equation  21. 
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(22)  L=£R(Pcg-Pl ) 


Pl  is  measured.  Pcg  is  estimated  using  the  measured  velocity  at  the  payload 
(Vl).  A  low-pass  filter  is  applied  to  Vl  to  remove  the  effect  of  the  oscillations.  Using 
this  filtered  velocity  (Vlf)  and  removing  the  wind  velocity  (Vw),  an  estimate  of  the 
velocity  at  the  center-of-gravity  is  obtained  (Vcg). 

(23)  VCG 

This  velocity  is  then  integrated  from  the  parachute  release  point.  The  estimated 
location  of  the  center-of-gravity  was  then  determined  using  Equation  22.  Figure  24 
presents  these  results. 


Figure  24.  Movement  of  Lever  Arm  Correction 


The  mean  lever  arm  corrections  are  Lx=0.0,  Ly=-0.1,  and  Lz=-3.94  feet.  Yavus 
and  Cockrell7  demonstrated  that  acceleration  of  the  air  mass  and  angle  of  attack  of  the 
parachute  significantly  affects  the  apparent  mass  coefficients.  Attempts  were  made  to 
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correlate  these  data  to  linear  and  angular  acceleration  of  the  payload.  However,  no  direct 
correlation  could  be  found.  Ideally,  acceleration  of  the  parachute  should  be  measured  in 
future  work.  With  the  existing  instrumentation  suite,  it  appears  the  flight  test  data  must 
be  filtered  to  remove  the  effects  of  oscillations  prior  to  parameter  estimation.  Using  the 
data  obtained  from  the  uncontrolled  test  drop  (3 1  ldropl),  a  filter  was  derived  using  the 
MATLAB®  System  Identification  Tools.  A  S^-order  low-pass  filter  with  a  cut-off 
frequency  of  0.009  (in  fractions  of  the  Nyquist  frequency)  was  selected.  Figure  25 
compares  the  measured  data  to  the  filtered  data  for  velocity. 


Figure  25.  Comparison  of  Measured  and  Filtered  Velocities 


This  filter  was  applied  to  the  position  and  velocity  data  for  all  drops  to  estimate 
the  motion  of  the  system  at  the  center-of-gravity.  The  obvious  concern  about  filtering  is 
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the  potential  of  eliminating  true  system  motion  while  removing  "noise."  Figure  26 
compares  the  glide  ratio  filtered  and  unfiltered  response  to  a  single  control  input. 


Implementation  of  the  designed  filter  appears  to  have  the  potential  of  eliminating 
desired  response  data.  The  character  of  the  response  (magnitude)  remains  in  the  filtered 
data  but  the  response  time  appears  to  be  adversely  affected.  A  filter  of  this  magnitude  is 
not  desired  but  may  be  necessary  for  parameter  estimation. 
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V. 


FLIGHT  TEST  RESULTS 


The  control  system  is  intended  to  affect  a  change  in  horizontal  velocity.  This  is 
best  demonstrated  by  assessing  the  glide  ratio  of  the  system  with  the  winds  removed. 
Figure  27  presents  the  glide  ratio  with  the  measured  control  inputs. 


The  results  show  that  a  nominal  glide  ratio  of  0.4  to  0.5  exists  with  no  control 
inputs.  Potential  causes  of  this  induced  glide  are  motion  induced  by  the  oscillations, 
imperfections  in  length  of  the  pressurized  actuators,  the  mathematics  of  creating  a 
horizontal  glide  ratio  which  eliminates  direction  of  motion,  or  errors  in  the  wind  estimate. 
This  nominal  glide  ratio  does  not  limit  the  assessment  of  the  response  due  to  control 
input,  as  we  are  interested  in  the  change  of  glide  ratio  at  the  time  of  control  activation.  At 
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time  zero,  all  PMAs  were  pressurized.  The  system  was  then  allowed  to  stabilize  to  a 
'trim'  condition.  A  change  in  glide  ratio  is  apparent  at  approximately  20  seconds  with  no 
change  in  the  state  of  the  controls.  The  first  incident  of  change  in  glide  ratio  can  be 
attributed  to  the  parachute  stabilization/inflation  process.  The  remaining  data  clearly 
show  a  correlation  of  glide  ratio  changes  to  the  activation  of  the  controls. 

A.  SINGLE  CONTROL  INPUT 

Figure  28  isolates  the  response  of  a  single  control  input.  An  increase  in  glide  ratio 
from  approximately  0.5  to  approximately  1.0  with  a  time  constant  of  about  4  to  5  seconds 
is  observed.  The  system  returns  to  its  oscillatory  trim  state  after  about  5  seconds 
following  removal  of  the  control  input.  The  reduced  magnitude  of  oscillation  or  coning 
angle  contributes  to  a  reduced  rate  of  descent  and  increased  glide  ratio. 


36 


B.  TWO  SIMULTANEOUS  CONTROL  INPUTS 


Recall  that  the  two  control  inputs  can  be  activated  simultaneously.  The  intent  is 
to  provide  additional  resolution  (every  45  degrees)  in  controlling  the  system.  Figure  29 
presents  the  response  with  two  simultaneous  inputs. 


As  exemplified  by  this  figure,  there  is  no  increase  in  performance  with  two 
control  inputs  over  that  achieved  with  one.  In  fact,  the  data  indicates  reduced  response 
results  from  the  simultaneous  activation  of  two  PMAs.  This  reduced  performance  is 
likely  due  to  leading  edge  collapse  (as  observed  in  the  ground-to-air  video,  Figure  30)  of 
the  parachute  with  two  control  inputs.  The  magnitude  of  the  oscillations  is  not  reduced  as 
dramatically  as  with  a  single  control  input. 
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Figure  30.  Parachute  Collapse  from  Two  Control  Inputs  (video  clip) 


The  vertical  velocity  of  the  system  was  observed  to  decrease  with  a  single  control 
input  but  remained  essentially  constant  with  two  simultaneous  control  inputs.  This 
reduction  is  likely  due  to  the  reduced  "coning"  angle  from  the  oscillation  damping  effect 
of  the  control  system.  Figure  3 1  presents  the  vertical  velocity  with  control  actuation  from 
a  single  and  a  double  input. 


Figure  31.  Vertical  Velocity  Response  to  Control  Input 


C.  ATTITUDE  RESPONSE 


For  control  system  design,  understanding  the  heading  and  yaw  rate  responses  of 
the  system  is  critical.  The  heading  and  attitude  rates,  as  measured  by  the  AHRS,  were 
correlated  to  the  control  inputs  for  analysis.  Actuation  of  the  control  system  resulted  in 
significantly  decreased  parachute  oscillations.  Figure  32  illustrates  the  roll  rate  data 
collected  from  the  AHRS.  Whenever  a  control  was  activated,  the  attitude  rates  were 
significantly  reduced. 


The  root  mean  square  (RMS)  of  the  roll  rate  reduced  from  16.8  degrees  per 
second  down  to  2. 1  degrees  per  second  with  control  activation.  This  behavior  was 
observed  each  time  a  control  was  activated. 
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A  plot  of  the  total  displacement  (or  coning)  angle  from  vertical  also  demonstrates 
the  reduction  in  oscillations. 


Figure  33.  Oscillation  Angle  with  Control  Input 
Figure  34  presents  heading,  measured  by  AHRS,  and  the  control  activation. 


Figure  34.  Heading  Response 
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With  careful  examination  of  these  data,  one  can  correlate  changes  in  heading  with 
control  inputs  as  indicated  with  the  arrows  on  the  graph.  However,  the  response  appears 
to  vary  greatly.  The  magnitude  and  direction  of  heading  changes  vary  with  the  activation 
of  the  same  control.  The  cause  of  this  variation  in  heading  response  has  not  yet  been 
identified.  Factors  such  as  the  coning  angle  at  the  time  of  control  input  must  be 
evaluated. 

During  initial  modeling  efforts,  it  was  apparent  that  the  rotation  rate  of  the  C-9 
parachute  severely  limited  the  efficiency  of  this  control  concept.  Excessive  control 
inputs  would  be  needed  to  affect  control  in  a  single  cardinal  direction.  Therefore,  it  was 
crucial  to  this  study  to  determine  if  the  operational  G-12  system  exhibited  similar 
characteristics.  The  results  from  this  demonstration  clearly  show  a  significantly  reduced 
rotation  in  the  larger  G-12  parachute  as  compared  to  the  C-9  parachute.  The  C-9 
rotations  can  be  observed  in  Figure  35.  . 


Figure  35.  C-9  Rotation  Rate  (video  clip) 
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Figure  36  presents  a  video  from  a  camera  installed  on  the  payload.  With  the  sun 
as  a  reference,  the  very  low  rotation  rate  of  this  system  can  be  observed.  A  ground  based 
video  (Figure  37)  also  demonstrates  a  low  rotation  rate  for  the  G-12  parachute. 


Figure  36.  G-12  Rotation  Rate  (video  clip) 


Figure  37.  G-12  Rotation  Rate  (video  clip) 
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VL  NAVIGATION  SENSOR  MODELS 


The  AGAS  is  expected  to  include  two  navigation  sensors:  1)  a  commercial  Global 
Positioning  System  (GPS)  receiver  for  position  determination  and  2)  a  heading  reference 
assumed  to  be  a  magnetic  compass  for  this  study.  To  assess  the  effects  of  navigation 
sensor  errors,  models  for  each  sensor  must  be  incorporated  into  the  simulation. 

A.  GLOBAL  POSITIONING  SYSTEM  (GPS) 

GPS  error  sources  include  errors  induced  by  the  atmosphere  (ionospheric  and 
tropospheric),  multi-path,  receiver  noise,  satellite  clock  noise,  and  Selective  Availability. 
Modeling  techniques  for  GPS  range  errors  resulting  from  these  sources  have  been 
developed  and  validated8  to  model  range  errors  and  not  errors  in  a  Cartesian  reference  as 
desired  here.  Cartesian  (x,  y,  z)  errors  would  therefore  have  to  be  formed  from  the  range 
errors  for  implementation  in  this  simulation.  This  necessitates  application  of  a  numerical 
solution  like  maximum  likelihood  techniques.9  Although  this  implementation  is 
relatively  trivial,  the  computation  resources  required  severely  limit  the  simulation  speed 
on  a  Personal  Computer  (PC).  Therefore,  a  variation  of  this  approach  was  implemented. 

1.  Selective  Availability 

Selective  Availability  is  a  means  of  intentionally  inducing  errors  into  the  GPS 
satellite  signal.  The  DoD  induces  these  errors  to  restrict  use  of  the  full  precision  of  GPS 
to  unauthorized  users.  Authorized  users  must  apply  a  receiver  capable  of  processing  the 
cryptographic  codes  to  remove  these  induced  errors.  Although  the  AGAS  concept  could 
incorporate  an  authorized  receiver,  it  is  desired  to  utilize  a  commercial  GPS  receiver  for 
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cost  savings.  Also,  with  airdrop  the  loads  may  not  be  fully  recoverable  and  loss  of  the 
authorized  receivers  would  not  be  desirable. 

The  errors  resulting  from  Selective  Availability  are  not  stochastic  in  nature. 
Therefore  system  identification  methods  were  employed  to  obtain  a  reasonable  error 
model.  Data  were  collected  at  the  Yuma  Proving  Ground  Satellite  Reference  Station.  An 
unauthorized  GPS  receiver  was  placed  on  a  known  survey  point.  GPS  position  data  were 
collected  for  over  two  hours.  The  position  data  in  the  three  Cartesian  axes  were 
differenced  with  the  surveyed  coordinates  resulting  in  the  Cartesian  errors.  These  errors 
represent  all  GPS  error  sources  identified  above.  Figure  38  illustrates  the  apparent 
random  nature  of  these  data. 
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Figure  38.  Measured  Unauthorized  GPS  Errors 

To  obtain  a  model  of  these  data,  the  MATLAB®  system  identification  toolbox 
was  utilized.  An  ARMAX10  model  was  utilized  with  the  input  being  white  noise  and  the 
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output  being  the  position  errors  shown  above.  The  ARMAX  model  incorporates  a 
prediction  error  method  with  a  model  represented  by  a  set  of  difference  equations  of  the 
form: 

(24)  A(q)y(t)  =  B(q)u(t  -  nk)  +  C ( q)e{t ) ;  where  y  and  u  are  the  outputs  and  inputs  of 
the  system,  respectively. 

The  coefficients  A,  B,  and  C  are  polynomials  that  describe  the  model's  difference 
equations.  The  prediction  error  is  minimized  using  an  iterative  Gauss-Newton  algorithm. 
The  ARMAX  function  returns  a  matrix  of  the  polynomial  coefficients.  This  matrix, 
referred  to  as  THETA  format,  can  then  be  transformed  into  a  transfer  function  using  the 
MATLAB®  command  TH2TF.  This  technique  resulted  in  the  following  transfer  function 
that  was  utilized  in  the  overall  system  model  to  obtain  GPS  errors: 

z4  -1.5302Z3  +0.2608z2  +  0.2566z  + 0.0192 
'  '  z4  -2.6500z3  +1.9582z2  +0.0337z  -0.3420 

Figure  39  presents  the  output  of  the  GPS  error  model  including  Selective 
Availability  Errors. 


Figure  39.  Modeled  GPS  (Unauthorized)  Errors 
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The  transfer  function  input  is  white  noise  initiated  with  a  random  seed  ensuring 
variable  errors  are  introduced  from  simulation  to  simulation.  To  assess  the  adequacy  of 
this  model,  the  mean,  standard  deviation,  and  root  mean  square  were  calculated  for  the 
measured  and  modeled  GPS  errors.  The  sample  of  measured  errors  presented  above  has 
a  mean  value  of  approximately  zero  meters  in  each  axis  and  a  standard  deviation  of  17.2, 
28.8,  and  21.1  meters  in  the  x-,  y-,  and  z-axes,  respectively.  The  modeled  results 
demonstrated  mean  errors  of  three  to  six  meters  with  standard  deviations  ranging  from  25 
to  35  meters.  The  root  mean  square  errors  for  the  model  were  found  to  range  from  26  to 
37  meters  for  three  independent  simulations.  The  model  produces  a  reasonable 
representation  of  the  measured  GPS  data. 

2.  Model  Without  Selective  Availability  Errors 

With  Selective  Availability  turned  off,  that  is  no  induced  errors,  a  commercial 
GPS  receiver  is  capable  of  navigating  with  greater  accuracy.  A  GPS  Error  Model  was 
derived  considering  a  noise  structure  proposed  by  Draper  Laboratory.11  This  report 
models  a  P-code  GPS  receiver  incorporated  into  the  Honeywell  Embedded  GPS/Inertial 
Navigation  System.  The  noise  model  incorporates  two  components,  accuracy  and  jitter. 
The  accuracy  noise  component  is  considered  exponentially  correlated  noise.  The  jitter 
component  consists  of  two  elements:  an  exponentially  correlated  noise  component  with  a 
faster  time  constant  than  the  accuracy  component  and  a  uniform  uncorrelated  noise 
component.  The  GPS  position  noise  model,  suggested  for  GPS-only  operations  (no 
inertial  aiding),  was  adapted  for  a  commercial  grade  (C/A-code)  receiver  by  adjusting  the 
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accuracy  and  jitter  standard  deviation  specification.  Table  3  presents  the  original  and 
adapted  models. 


Draper  P-code 
Model 

C/A-code 

Model 

Accuracy  Standard  Deviation 

30 

45 

Accuracy  Correlation  Time  Constant  (seconds) 

0.1 

0.1 

Jitter  Standard  Deviation  (feet) 

5 

8 

Jitter  Correlation  Time  Constant  (seconds) 

0.05 

0.05 

Uniform  Uncorrelated  Noise  Standard  Deviation 

q/(2*sqrt(3)) 

q/(2*sqrt(3)) 

Note:  q  is  the  quantization  interval 


Table  3 .  GPS  Position  Error  Model 


This  model  was  incorporated  into  the  Simulink®  simulation.  Figures  40  and  41 
illustrate  the  results  obtained  from  this  model  for  a  400-second  simulation  and  50-second 
simulation,  respectively.  The  standard  deviation  of  the  3-axis  error  for  this  simulation 
was  56.3  feet  (17.2  meters)  which  is  close  to  the  specifications  for  an  commercial 
receiver  with  Selective  Availability  off. 


Figure  40.  GPS  Error  Model  -  SA  Off 
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Figure  41.  GPS  Error  Model  -  SA  Off  (50-second  Simulation) 


B.  HEADING  SENSOR 

The  heading  sensor  is  assumed  to  be  a  magnetic  compass  for  this  study.  Two 
components  of  errors  are  considered  here,  a  static  error  or  bias  and  a  dynamic  (noise) 
component.  System  specifications  for  the  Attitude  Heading  Reference  System  (AHRS) 
provide  a  static  error  of  +2  degrees  (+  1  degree  with  velocity  aiding)  and  a  dynamic 
component  of  +2  percent.  The  AHRS  incorporates  rate-gyros  to  obtain  3-axis  attitude 
rates  and  attitude  data.  Specification  sheets  of  a  low-cost  digital  magnetic  compass 
produced  by  KVH  Industries  presented  similar  accuracy  statements.  The  static  error  is 
incorporated  as  a  bias  element  in  the  Simulink®  model  and  is  set  as  a  uniform  random 
variable  at  the  start  of  each  simulation.  The  dynamic  component  is  found  by  adding  2% 
of  the  current  heading  reading.  Figure  42  presents  a  400-second  simulation  of  the 
heading  error. 
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Vn.  EQUATIONS  OF  MOTION 


The  kinematic  equations  of  motion  are  determined  assuming  the  parachute 
and  payload  system  are  a  rigid  body.  A  six-degree  of  freedom  model  is  considered  with 
the  motion  about  the  center  of  gravity.  The  notation  and  derivation  follow  the  convention 
presented  by  Isaac  Kaminer12  as  well  as  by  Antonio  Pascoal  and  Carlos  Silvestre13.  The 
derivation  first  considers  an  origin  displaced  from  the  system's  center  of  gravity  by  a 
distance  P0.  In  the  implementation  of  these  equations  to  this  problem,  the  origin  is 
moved  to  the  center  of  gravity  of  the  system. 

A.  NOTATION 

The  following  illustrates  the  reference  system  used  for  this  effort: 


Figure  43.  Coordinate  Convention 
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The  following  notations  are  used  in  the  development  of  the  equations  of 

motion. 

-  Capital  letters/symbols  denote  vectors  or  matrices 

-  Small  letters/symbols  denote  scalars 

-  {B}  \  name  of  the  coordinate  system,  B  denotes  body-axis,  W  denotes  wind- 
axis,  and  U  denotes  inertial  axis. 

-  caR  :  rotation  matrix  from  coordinate  system  {A}  to  coordinate  system  {C} 

-  CPD :  position  of  point  D,  measured  in  {C}  and  expressed  in  {C} 

"  A(C  Pd)  '■  position  of  point  D,  measured  in  {C}  and  expressed  in  {A}  where 

a(cpd)=ac  R(cpd)*cpd 

-  CVD :  velocity  of  point  D,  measured  in  {C}  and  expressed  in  {C} 

-  cQd:  angular  velocity  of  point  D,  measured  in  {C}  and  expressed  in  {C} 
d 

-  — :  time  derivatives  in  the  body-axis  {B} 

-  (' ):  time  derivatives  in  the  inertial-axis  {U} 

NOTE:  if  the  symbol  for  the  coordinate  axis  is  omitted,  the  inertial-axis 
{ V }  is  assumed. 

B.  ASSUMPTIONS 

Throughout  this  effort,  the  following  assumptions  were  considered: 

-  Rigid-body  system 

-  Only  the  motion  after  complete  parachute  deployment  is  considered. 

-  Non-rotating  earth,  i.e.  ECEF  coordinate  equals  inertial  coordinate 

-  Wind  axis  equals  inertial  axis 

C.  DERIVATION  OF  EQUATIONS  OF  MOTION 

Starting  with  Newton's  Law,  F  =  ma 
uF=muVB 

The  inertial  velocity  of  the  body  is: 

UVB  =UbP{BPb  +^x^0}  ;  where  P0  is  the  position  of  the  parachute 

reference  to  the  origin  (center-of-gravity). 
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Taking  the  derivative: 

(26)  UVb=ubr\^bVb  +QxP0  +QxPq ^+vbr(bVb  +  fix PQ ) ;  where  PQ  =  0  due  to 
rigid  body  assumption 

recall: 

~P 

(27)  ubR=ubR-S(b nB),  bQ.b  =  q 

r 

0  -r  q 

and  S(bQ.b)=  r  0  -p  (skew  symmetric  matrix) 

-q  P  0_ 

Substituting  (27)  into  (26): 

(28)  uvB="RpL‘v,  +  n  x  />0|+"i!  ■  s(‘n,  +  n  x  ) 

Noting: 

(29)  s(BnBfvB=BnBxBvB 

Substituting  (29)  into  (28): 

vvB  =%R&r,  +  n  x  P0  |+^s  n,x(evt  +cixP0) 


Rearranging  Terms: 

-P0  x  xfi+Y.) 

Consider  forces  in  the  body  axis: 

bF=bR-uF=bR-MtuVb;  where:  Mr  =m/3 
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( 33)  U  Lb  =uNb  - Mt  Pq  xu  Vb  ;  angular  momentum  equals  total  applied  moment; 

where:  M  =  ml.  ;  I3  is  a  3x3  identity  matrix  and  m  is  the  system's  mass 
Recall: 


uvB=uBRBvB 

"V^'iR^e+lk’V, 

iR’V^Blpx’V,) 

(34)  vV,=lR^BVB+asIl(axBV,) 

Substituting  (34)  into  (33): 

ulb=vn„  -mtp„  j 
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(36) 


(37) 


(38) 


uLb=ubRbL ;  nB=lRjBL+uBRBL=lR^BL+uBR{QxBL) 
(35)  ubR^bL+ubr(Qxbl)=uNb  -MTP0  xf  ubR^bVb+ubr(QxbVb)  ) 


dt 

By  definition: 


3L  =  IbQ  ;  where  Is  is  the  system's  moment  of  inertia  matrix 


BL  =  IBQ 


Substituting  (36)  into  (35) 


uBRiBn+uBR{nxiBn)=uNB  -mtp0  x  ubr^-bvb+ubr(qxbvb) 

dt 


Rearranging  Terms: 


.  uNb  =vBRlBQ.+uBR(a X IbQ) + MtP0  x^RjtBVB+uBR.(QxBVB) j 
Transforming  into  the  body  axis  {B}: 

bn,=;runb  =i,  3- n+(nxiBn)+MTp0  xfl'r,  +(nx*r,)J 


"nb  =  iB)-n+(nxiBn)+MTp0  x(n><‘vB)+MTp0  x{^-‘vB 

dt  \dt 

Applying  (31)  to  (37): 


bNb  =Mr—Q+M 
dt 

Rearranging: 


RT 


f  rf  \ 

—  vB 

\dt  Bj 


■  Q  x  MRn + Mrt  (fix5  Vb  ) 


j + n  xMfi + (fix*  V,  )| +mb-'  bn. 


Combining  (38)  and  (32)  in  state  space  form: 
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Mt-t%  +MmQ+BnB  X  \Mmn  +Mt  bVb  } 


sj  }+QxMRn+MRT{nxBvB) 


dr, 
*V‘ 
— n 

.  dt 


^r~!k*6+*ns  x^q+a/,' bvb}} 


Mr~)Mzt  ~t*Vb  +&xMRn+MRT{QxBVB) 


0  mSI’n, 


D.  APPARENT  MASS 

As  a  body  accelerates  through  a  fluid,  the  fluid  itself  must  accelerate  to 
accommodate  the  motion  of  the  body.  Resultant  forces  and  moments  are  applied  to  the 
body.  A. common  method  for  accounting  for  these  forces  is  to  include  added  or  apparent 
mass  terms  in  the  equations  of  motion.  Sir  Horace  Lamb  performed  the  original  work  on 
the  effects  of  accelerating  fluid  on  a  body. 14  Lamb  derives  the  apparent  mass  effects  and 
identifies  15  independent  terms.  Cockrell  and  Doherr15  showed  that  for  a  revolution  of 
the  body  about  a  plane  of  symmetry,  the  number  of  independent  terms  reduces  to  the 
following: 

ocll  -  motion  along  the  x-axis 
a -  motion  along  the  y-axis 
a33  -  motion  along  the  z-axis 
#44  -  motion  about  the  x-axis 
a55  -  motion  about  the  y-axis 

a\ s  ‘  coupling  between  motion  along  the  x-axis  and  about  the  y-axis. 
a24  -  coupling  between  motion  along  the  y-axis  and  about  the  x-axis 

These  terms  have  been  demonstrated  to  change  significantly  with 
acceleration  of  the  body.16  However,  estimation  of  these  changes  is  very  difficult  and 
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beyond  the  scope  of  this  research.  Therefore,  only  the  scalar  values  of  apparent  mass  and 
moment  of  inertia  terms  are  considered. 

1.  Estimation  of  Apparent  Mass  Terms 

Doherr  and  Saliaris17  demonstrated  that  for  an  equivalent  parachute,  four 
independent  apparent  mass  terms  exist: 

A  D  - 

(40)  <*33  =  {—  )3}  ;  Dp  is  the  profile  diameter  of  the  parachute 

3  2 

(41)  an  =  #22  =  2^33 

(42)  ccl5  =  —a.  24  =  0.2ojj  fb 

(43)  <*44  =  a55  =  0.048 Dp2au 

2.  Effects  on  the  Equations  of  Motion 


The  force  imparted  on  the  body  by  the  air  mass  can  be  determined  by  starting 
with  the  equations  of  motion  for  the  body  itself  (Equation  39)  and  replacing  the  mass 
matrices  with  apparent  mass  matrices  as  shown  here: 


Mt  Af,  Mtr  Ajr,  Mrt  "9  -Am',  Mr  ^  Ar 


(44) 


(45) 


Mt  —bVb  +MtrC1+b Qb  X  {mtrq +mt  bvb  } 
dt 

M,3rn+Mj(-‘vB\nxMKn+M„{n*’v,) 

at  \Cu  J 


<*„  0  0  <*44  0  0 

AT  =  0  a 22  0  ;  Ar  =  0  <*55  0 

0  0  <*33  0  0  0 
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where: 

AdT  Ad ' jR  Mr  "T"  Aq *  AvI ~f~  Aqg 

AA Rt  A/I R  ~  Ad RT  "*■  ■Art  Mr  ^r 
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E.  EXTERNAL  FORCES  AND  MOMENTS 


The  external  forces  and  moments  on  the  body  include  the  effects  of  aerodynamics 
(including  actuation  of  the  control  surfaces)  and  gravity. 

=  FaERO  +  Fm  +  FgRAV 

Fb  =  FfAERO  +  Nm  +  Ngrav 

Buoyancy  terms  as  discussed  by  Gockel18  are  not  included  as  the  convention  used 
herein  does  not  include  the  mass  of  the  trapped  air  in  the  overall  systems  mass.  As 
Gockel  demonstrates,  if  one  is  to  include  the  mass  of  the  trapped  air  a  compensating 
"buoyancy"  term  must  be  added  to  the  equations  of  motion. 

1.  Aerodynamic  Forces  and  Moments 

Using  a  first  order  Taylor  series  expansion  about  a  nominal  trim  condition  and 
following  Kaminer's19  notation,  the  aerodynamic  effects  are  found  to  be: 

F aero  -Fq+  SFx  +  SFxx  +  <5Fa A 

Naero  ~  N0  +5NX  +SNix  +  SNAA 

where  x  and  x  are  the  state  vector  and  first  derivative  of  the  state  vector  and  A  is 
the  control  vector.  5F  and  5M  are  the  partial  derivatives  of  the  forces  and  moments  with 
respect  to  each  parameter  in  the  state  vector  and  their  first  time  derivative. 

Non-dimensionalizing  and  putting  in  state  space  form: 

F aero  ~  $JfcFo  +FCFx  +SCFix+SCFAA] 

F AERO  —  qS0dQ  {Cyo  +SCNx  +SCNix  +  SCNaA} 
where: 
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cos^cos#  sin^cos#  -sin# 

where  br_  C0Sy/  sjn$  sin ^  -  sin cos ^  sin^sin^sin^  +  cos^cos^  cos^sin^ 
cos^sintfcos^+siny/’sin#  sin^sintfcos^ -cosy/- sin  ^  cos#  cos^ 

Recall  that  P0  is  defined  as  the  distance  from  the  origin  to  the  center  of  gravity. 
The  gravitational  force  creates  a  moment  about  the  defined  origin. 

0  " 

0  ;  later  P0  will  be  defined. 

mg 

3.  Total  External  Forces  and  Moments 

The  total  external  forces  and  moments  are  now  given  as: 

'  0 

BFB=qS,\SCFo+SCp,+SCrix+SCF,A\+;R  0 

mg 

1  r°" 

mg 

In  matrix  form: 


bFb  _  qS0I3  0  {5CFo+SCFx+SCFix+SCFiiA}  +  BR 
_B N B _  _  0  qS0daI3_  J>ou^_ 

F.  COMPLETE  EQUATIONS  OF  MOTION 

This  derivation  assumes  that  the  system  is  a  single  rigid  body  and  resulted  in  a 
general  set  of  equations  of  motion.  These  equations  could  then  be  applied  to  the 
parachute  and  payload  separately.  For  example,  to  fully  capture  the  oscillations  of  the 


61 


parachute  system,  it  may  be  beneficial  to  establish  a  coordinate  system  at  the  parachute 
aerodynamic  center.  The  model  would  then  incorporate  the  applicable  components  of 
these  equations  for  both  the  parachute  and  payload.  The  general  equations  of  motion  are 
completely  defined  as: 

£‘r,  1  r  K~'lV'T,h+,n,xfa'mn+M'‘vB}} 

+  \K~'  o  lr?5„/3  0  J&CFl,+SCrx+6CFii  +  SCFiA}']  r  ‘R ' 

L  o  «r'JL  ° 

where: 


~  M’r 

M'm 

iSArj,  4*  Ay 

M  tr  4  Atr  1 

Mrt 

Mr. 

Mrj  4  Art 

4 

i _ 
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Vm.  ANALYSIS  OF  EQUATIONS  OF  MOTIONS 


A.  EQUATIONS  OF  MOTION  EXPANSION 

Expansion  of  the  equations  of  motion  is  necessary  for  parameter  estimation.  This 
expansion  also  serves  to  validate  the  derivation  of  these  equations  against  existing  works 
in  the  field.  First,  the  coordinate  system  axis  is  selected  aligned  with  the  parachute  body 
axis  but  displaced  from  the  center  of  gravity  by  a  distance  z.  Expanding  the  equations  of 
motion: 

X  =  (m  +  an)u  + ( mz + ax  5  )q + (m  -  a24  )pr  -  (m + a22  )rv + (m + a33  )qw 

Y  =  (m+a22)v-  (mz  -a24)p  +  (m+an)ru  +  (mz +al5)qr-(m+a~3)pw 
Z  =  (m  +  a~3)w  -  (mz  +  al5)q2  -  (m  +  an)qu -(m-a2A)p2  +(w  +  a22)/?v 

L  =  (IX  +aM)p-  (mz  -  «24  )v  +  (1 2  -Iy-  ass)qr  -  (mz  +  al5)ru  -  (mz  -  a24)pw  +  (au  -cc~)uw 
M  =  (Iy  +a5S)q+(mz+al5)u  +  (Ix -I2  +aA4)pr  ~(mz+a,5)(qw-rv)-(an  -a33)uw 
N  =  /Zr 

These  equations  are  validated  by  the  work  of  Cockrell  and  Doherr20. 

Taking  terms  like  rv,  qw,  etc.  to  be  small  and  eliminating  non-linear  terms  like  p2: 

X  =  (m+au)u  +  (mz+cc\5)q 

Y  =  (m  +  a22)v  -  (mz  -  a24)p  . 

Z  =  (m+a33)w 

L  =  (/*  +a44)p-(mz -a24)v 
M  =  (Iy+  oc55)q  +  (mz  +  al5)u 
N  =  1/ 

Forming  forces  and  moments  as  an  initial  (trim)  value  plus  a  perturbation  and 
applying  gravitational  forces  and  moments. 

X  =  AX-mgsin(0)  =  (m  +  auyu  +  (mz+ccl5)q 

Y  =  XY -mgcos(0)sin(^)  =  (m+a22)v-(mz-a2i)p 
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Z  =  Z0  +  AZ  +  mg  cos (0)  cos(0)  =  (m  +  a33)w 
L  =  AL  +  mgz  sin(^)  =  (7X  -sraAA)p  -  ( mz  -  a24)v 
M  =  AM  -  mgz  sin(0)  =  (/^  +  a55  )#+ (mz  +  a,  5  )w 
N  =  AN  =  I/ 

Subtracting  initial  conditions  as  defined  as  X0=0,  Z0=-mg,  Y0=0,  L0=0,  Mo=0, 
N0=0,  0o=O,  <J)0=0  and  applying  small  angle  assumption  (e.g.  cos  6  «l,sin#«6>)  two  sets 
(Longitudinal  and  Lateral-Directional)  of  uncoupled  equations  are  formed. 

B.  LONGITUDINAL  EQUATIONS 

(m+an)u  +  (mz+aX5)q  =  AX-mgO 
(m+a33)w  =  AZ 

(mz+aX5)u  +  (Iy+a55)q  =  AM-mgz6 


The  external  forces  and  moments  are  approximated  with  a  Taylor-Series 
expansion: 

AX^AX.^  +  AXt^+AX^  +  AX^  +  AXj  +  AXj  +  AXtS 


Ar=AZ*F+A2‘^-+AZ--^+AZ*^+A^+AZw+AZ^ 


wn 


wn 


AM  -  iT  +  AM*  w  +  +  AM*  +  +  &Mss 


W0  UW0  ”W0 


Reducing  terms  that  have  no  apparent  effect  on  forces  or  moments  in  a  particular 
direction  (e.g.  AA^w,  A Z-u)  and  substituting: 


W0(m+au)^  +  (mz+al5)q  =  AXu^r  +  AX.^r-mgO+AXsS 

W  o  Wq  yy  0 


«/  1J  JJ  W  w 

rr0(/KZ +a15)—  +  (/,  +a5S)$  =  AM,  —  +  AM,  —  +  AA/„  — +  AMW  — +  A Mqq  -  /Kgz0  +  A MsS 


Defining  the  stability  derivatives: 

AX.  „  AZ; 


X  = 


;  z,=- 


W0{m+any  '  Pf0(/M+a33) 

where:  i  =  u,u,w,w 


;  M( 


AM. 


*W+«55> 


x.= 


AX,. 


AZ, 


AMy 


' _ •  z.= _ J- _ •  M .  =  • 

'  0  +  «n)’  '  (m  +  a33)’  ;  {Iy+a  55) 

where:  j-q,q,8 

Establishing  the  coordinate  system  now  at  the  center  of  mass  of  the  system, 
substituting  the  stability  derivatives  and  rearranging  terms: 

Xx8 


V  u  '  ITT 


a , 


15 


W0(m+au) 


n-Y  U  -  m%9 

$  TIr  ] 


W0  W0(m+an)  W0 


(i-z  )—  =  z  ZL+^A+M. 

1  *V„  -Wc  W0  W0 


w„ 


«1 


15 


-M, 


J 


K(Iy+cc55) 

In  state  space  form  with  0  =  q: 


—  -  WM.  —  +  q  =  WqM  —  +  WMw  —+Mq+Ms8 
W0  °  w  Wo  Wo  *  Wo  q  S 


w, 


(l-x4) 

0 

ais 

(I  y  +a55) 
0 


a,, 


0  — - t  0 

W0{m+au) 

(1-Zj  0  0 


-Mi  |  -WM* 
0 


u 

w 

w. 

4 

[  #  _ 

0 


0  ~mg  1 

U 

pq 

W,{m+an) 

K 

K 

—  0 

w 

+ 

z* 

K 

Wo 

K 

M f  0 

9 

Ms 

1  0 

6 

0 

C.  LATERAL-DIRECTIONAL  EQUATIONS 

{m  +  a^u  -  (mz  -a2A)p  =  AY-  mg<j> 

-  (i mz -a2A)v  +  (Ix  +au)p  =  AL  + mgz<j) 

I/  =  AN 


65 


The  external  forces  and  moments  are  approximated  with  a  Taylor-Series 


expansion: 


AF  =  A7V  —  +  A  Y,  —  +  A  YpP  +  A Ypp  +  A  Yrr  +  AY.r  +  A  YsS 


&L  =  ALV  +  AL,  +  ALpp  +  ALpp  +  ALrr  +  AL.r  +  ALsS 


AA/-  =  ANV  +  AN+  + ANpp  +  ANpp  +  ANrr  +  AN-r  +  ANsS 

W  o  W o 

Reducing  terms  that  have  no  apparent  effect  on  forces  or  moments  in  a  particular 
direction  (e.g.  AYpp,AN,v)  and  substituting: 

W0(m  +  a22)^r-  (mz  -  a2A)p  -  AYV  ~  +  AY,  +  A Yrr  +  A YsS  -  mg<f> 

yy  0  yv0  w 0 

~Wg{mz  -  au)  —  +  {I x  +a44)p  =  ALV  +  AL,  +  ALpp  +  ALrr  +  ALS5  +  mgz<j> 


W0  pl 


I2r  =  ANV  —  +  ANpp  +  ANrr  +  ANsS 
"o 


Defining  the  stability  derivatives: 


y  _  AY.t 
W0(m+a22)  ’ 


wc(L  +«„)  ’  ’  vj. 


where:  /  =  v.v 


L.=-^l  .  ■ 

*  ST  .  _  \  »  r 


(tn + cc22)  {Ix+cc^) 


where:  j  =  p,p,r,r,S 

Establishing  the  coordinate  system  now  at  the  center  of  mass  of  the  system, 
substituting  the  stability  derivatives  and  rearranging  terms: 
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D.  SUMMARY  OF  EQUATIONS 

In  summary,  two  sets  of  independent  equations  have  been  derived  for  the 
longitudinal  and  the  lateral-directional  motion.  These  equations  form  the  models  to  be 
used  for  parameter  estimation. 

Longitudinal: 


Lateral-Directional : 


(l-Tj 

#24 

Ix  +CC  44 
0 

0 
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Non-dimensionalizing : 

F  =  qS0CF  and  M  =  qS0d0CM ;  where:  q ,  S0,  and  do  are  the  dynamic  pressure, 
parachute  reference  area,  and  the  parachute  reference  diameter,  respectively. 
Longitudinal: 


Which  are  in  the  form  Inx  =  A„x + Bnu ;  define:  A  =  In~xAn  and  B  =  In~xBn  to  get 
in  the  form:  x  =  Ax +  Bu 
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IX.  PARAMETER  ESTIMATION 


A.  OVERVIEW 

Parameter  estimation  is  the  process  of  determining  aerodynamic  stability  and 
control  derivatives  from  flight  test  data.  A  large  number  of  techniques  have  been 
employed  for  parameter  estimation  in  aircraft  research.  This  study  investigated  the 
application  of  Maximum  Likelihood  Estimation  methods  described  by  Richard  Maine 
and  Kenneth  Iliff  of  NASA21. 

Wind  tunnel  data  on  the  stability  and  control  derivatives  for  this  system  are  not 
available.  This  effort  encompassed  an  initial  start  at  tackling  this  problem.  Parachute 
aerodynamics  are  truly  non-linear  and  significantly  coupled  across  axes.  The  maximum 
likelihood  techniques  implement  linear  models  with  no  cross-coupling. 

B.  MODIFIED  MAXIMUM  LIKELIHOOD  ESTIMATION 

Once  the  system  is  modeled,  the  system's  equations  of  motion,  expressed  in  state 
space  form,  are  utilized  to  model  the  dynamic  response  due  to  a  simulated  input  with  an 
initial  "guess"  of  the  aerodynamic  parameters.  The  modeled  output  is  compared  to  the 
flight  test  results;  the  unknown  parameters  are  adjusted  as  to  minimize  the  output  error. 
Detailed  descriptions  of  the  MMLE3  algorithm  as  implemented  in  the  MATLAB®  are 
provided  in  the  user's  guide22.  A  block  diagram  of  this  concept  is  provided  in  Figure  44. 

A  coning  motion  characterizes  the  dynamics  of  a  flat  circular  parachute.  This 
motion  will  likely  not  be  identified  during  the  estimation  process  due  to  the  averaging 


®  MATLAB  is  a  registered  trademark  of  Mathworks,  Inc.,  Natick,  MA 
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techniques  implemented.  Karl  Doherr  suggests23  the  addition  of  noise  to  the  average 
normal  force  or  pitching  moment  coefficients  to  compensate  for  this. 


Figure  44.  Parameter  Estimation  Concept 


C.  RESULTS 

The  MMLE3  algorithm  as  implemented  in  MATLAB386  (an  obsolete  version  of 
MATLAB®)  was  successfully  run  for  several  test  cases.  First,  MATLAJB386  was  hosted 
on  a  Windows  95®  on  a  PC.  The  program  needed  to  be  run  in  the  DOS  mode  and  could 
not  be  run  in  a  Windows  shell.  In  addition,  several  memory  management  conflicts  had  to 
be  resolved  before  being  able  to  run  the  program.  After  successfully  launching  the 
program,  the  Naval  Postgraduate  School's  (NPS)  MMLE3  implementation  for  parameter 
estimation  was  executed.  Several  test  cases  were  executed  including  the  T37  and  F14 
aircraft  simulations.  In  all  cases,  duplicate  results  to  those  obtained  by  Graham24  were 
obtained. 


®  Registered  Trademark  of  the  Microsoft  Corporation 


70 


The  NPS  implementation  was  then  modified  to  incorporate  the  equations  of 
motion  for  the  parachute  (identified  above).  Simulated  data  (using  estimates  for  the 
stability  derivatives)  were  generated  and  qualitatively  observed  to  have  the  same  general 
character  of  the  flight  test  data.  However  when  utilizing  the  simulated  data  to  estimate 
the  applicable  stability  derivatives,  the  algorithms  failed  to  predict  the  stability 
derivatives  utilized  to  generate  the  simulated  data.  The  program's  execution  would  halt 
after  receiving  errors  indicating  the  HES  matrix  was  near  singular.  This  matrix  rotates 
and  scales  the  gradient  to  provide  a  single  step  convergence  for  the  quadratic  cost 
function.25  Initially,  it  was  believed  that  this  error  resulted  from  the  model  being 
overparameterized.  After  significant  investigation  and  numerous  trials  using  simulated 
and  flight-test  data,  no  results  could  be  obtained  using  these  algorithms. 

The  implementation  of  the  F-14  models  was  then  analyzed  again.  The  F-14 
equations  of  motion  were  utilized  and  the  F-14  parameters  were  modified  to  reflect  a 
slower  aircraft  (approximately  the  speed  of  the  descent  rate  of  the  parachute).  Simulated 
data  were  successfully  generated  but  when  applying  the  MMLE3  algorithm,  the  HES 
matrix  again  went  singular.  It  was  concluded  that  the  MMLE3  algorithm,  as 
implemented  in  MATLAB386,  could  not  provide  estimates  of  parachute  stability 
derivatives.  There  appears  to  be  a  numerical  sensitivity  in  the  algorithms  to  the  relatively 
low  numbers  representing  the  parachute  system.  Further  work  may  result  in  the 
successful  application  of  this  algorithm  to  the  parachute  problem.  However,  this  is 
beyond  the  scope  of  this  study.  It  is  also  recognized  that  the  utility  of  results,  should  they 
be  achieved,  may  not  be  sufficient.  This  implementation  requires  linearization  of  the 
equations  of  motion.  The  parachute  dynamic  behavior  is  clearly  not  linear.  In  addition. 
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other  research26  has  demonstrated  that  the  stability  derivatives  themselves  are  functions 
of  angle  of  attack,  and  may  even  vary  with  time.  Figure  45  shows  how  the  drag,  lift,  and 
moment  coefficients  change  for  a  typical  flat-circular  parachute. 


Figure  45.  Flat  Circular  Aerodynamic  Coefficients 
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X.  SYSTEM  MODEL 


Without  successful  MMLE  results  and  in  the  absence  of  sufficient  wind  tunnel 
data  for  the  C-9  parachute,  a  point-mass  system  was  assumed  with  the  only  forces  on  the 
system  being  drag  and  weight.  The  applicable  equations  of  motion  are: 

(50)  X-  (m  +  an)u  = -D  cosy  cosiff 

(51)  Y  =  (m  +  a22)v  =  -D  cos y  sin  if/ 

(52)  Z  =  (m  +  a~3)w=-Dsmy  +  W 

where:  W  is  the  calibration  system  weight,  y  and  if/  are  the  flight  path  angle  and 
yaw  angle  respectfully.  D  =  qCDS  D  is  drag. 


(53) 

(54) 

(55) 

mass  (airspeed),  Vw=wind  velocity.  Assumes  no  rotation  between  the  fixed  earth 
reference  and  the  system's  body  axis. 

VG=VA+K 


The  reference  angles  are 


w  Jv2  -  w2  .  v  u 

sm/  =  — ;  cos y  =  - — - ;  sin^=  —  =  -=— ;  cos y/  =  - 


VT 


■yJVf  ~  W2 


Substituting,  rearranging  terms,  and  putting  in  state  space  form: 
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where:  are  the  apparent  mass  terms,  here  assumed  to  he  constant. 

Noting  that: 

VG  -VA+VW\  where  Vg  =  ground  velocity,  VA=velocity  relative  to  the  local  air 
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Equation  (54)  presents  the  equation  for  VA .  Substituting,  the  EOM  for  estimating 
ground  speed  results. 

w+au  o  o  t1  r  r«i  r  o  T 

0  m  +  ccn  0  <  ~-q(y°S  v  +  0  >  +  Vw 

0  0  m  +  a3.\  [  T  [_wj  [w\ 

Equation  (56)  can  be  solved  numerically  to  estimate  the  system  response.  To 
estimate  heading,  a  constant  rate  of  rotation  was  assumed.  The  flight  test  results  for  the 
C-9  parachute  showed  significant  rotations  and  changes  in  rotation  rate  with  control 
activation.  Flight  testing  of  the  larger  G-12  parachute  showed  a  mean  rotation  rate  of 
approximately  1.8  degrees  per  second.  A  standard  deviation  of  1.0  degree  per  second 
was  assumed  based  on  qualitative  observations  during  flight  test.  An  insufficient  sample 
size  was  available  to  accurately  determine  a  standard  deviation  for  rotation  rate.  A 
normal  random  number  generator  at  the  start  of  each  simulation  determines  the  rotation 
rate.  That  rate  is  then  integrated  to  provide  heading.  This  model  was  incorporated  into 
Simulink®.  Multiple  wind  profiles  were  incorporated  into  the  simulation  using  a  lookup 
table  based  on  the  current  system  altitude.  The  equations  of  motion  are  contained  in  a 
MATLAB®  function  block  which  calls  a  script  file  titled  DOF3.  This  script  file 
(Appendix  E)  also  incorporates  the  parachute  system  parameters  and  apparent  mass 
terms. 

Figure  47  presents  the  measured  velocity  data  from  flight  test  as  compared  to  the 
modeled  velocity  data  for  an  uncontrolled  drop.  The  'noise'  in  the  measured  data  results 
from  the  velocity  being  measured  at  the  payload  which  is  experiencing  significant 
oscillations.  Since  the  point  mass  model  does  not  incorporate  these  oscillations,  no 
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'noise'  is  apparent  in  the  modeled  data.  As  is  demonstrated  in  the  graph,  the  velocity  data 
agree  very  well  for  this  uncontrolled  condition.  For  this  run,  atmospheric  density 
measurements  were  not  available.  Therefore,  the  modeled  descent  rate  does  not  match 
precisely  with  the  measured  descent  rate  but  the  differences  appear  negligible. 


Figure  46.  Simulink  Realization  of  System  Model 


Figure  47.  Measured  vs.  Modeled  Velocity 


Initializing  the  model  at  the  start  position  of  the  flight  test,  the  model's  ability  to 
estimate  position  of  the  system  was  evaluated.  Figure  48  shows  the  model  very 
accurately  predicts  the  flight  path  of  the  system  under  the  given  wind  conditions. 


Figure  48.  Measured  vs.  Modeled  Position 


The  response  of  the  system  to  a  control  input  was  identified  during  flight  test.  For 
a  single  control  input,  a  glide  ratio  of  0.4  to  0.5  was  found  with  a  time  constant  of  about 
4-5  seconds.  The  glide  ratio  was  decreased  to  approximately  0.2  when  two  simultaneous 
controls  were  actuated.  This  response  was  incorporated  into  the  system  model  by  adding 
the  resultant  horizontal  velocity  components.  Recall  the  flight  test  results  demonstrated 
an  apparent  decrease  in  descent  rate  upon  control  activation  resulting  from  the  decreased 
oscillation  angles.  This  reduced  vertical  velocity  is  an  apparent  benefit  to  the  system 
allowing  more  time  in  the  air  to  be  controlled  to  the  desired  impact  point.  However,  the 
implementation  of  the  point  mass  model  does  not  predict  oscillation  angle.  Therefore,  no 
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vertical  velocity  change  was  incorporated  into  the  model  on  control  input.  The  modeled 
response  to  double  and  single  control  inputs  is  illustrated  in  Figure  49. 


1  control 
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Figure  49.  Modeled  Control  Input  Response 
As  observed  in  flight  test  (Figure  50),  a  glide  ratio  of  approximately  0.2  for  two 
simultaneous  control  inputs  and  0.4  to  0.5  for  a  single  control  input  result.  A  time 
constant  of  approximately  3  to  5  seconds  is  observed. 


Figure  50.  Measured  Control  Input  Response  (Two  Controls  and  Single  Control) 


The  modeled  data  closely  represents  the  measured  system  response. 
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XL  SIMULATION 


A.  OVERVIEW 

The  individual  modeling  efforts  have  been  presented.  Three  major  components 
are  included  in  the  overall  system  model:  1)  Dynamics  Model,  2)  Sensor  Model,  and  3) 
Control  System  Model.  In  addition,  the  reference  trajectory  generator  was  implemented 
using  the  same  equations  of  motion  utilized  in  the  dynamics  model.  This  trajectory 
generator  has  the  additional  benefit  of  being  a  tool  for  calculating  the  Computed  Air 
Release  Point.  Figure  51  provides  a  block  diagram  of  the  overall  control  concept.  This 
control  methodology  was  implemented  in  Simulink®  (Appendix  C)  and  integrated  with 
the  sensor  models  and  dynamics  model.  Single  simulation  runs  can  be  executed.  It  is 
desirable  to  assess  to  obtain  a  statistical  base  for  many  trials. 


Figure  51.  Control  Concept 
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B.  INITIALIZATION 


To  obtain  the  statistical  base  desired,  an  initialization  file  was  created  where 
numerous  parameters  utilized  in  the  simulation  were  selected  randomly.  The 
initialization  file  then  provided  the  parameters  required  for  numerous  simulations.  A 
MATLAB®  script  file  called  INIT  was  created  to  select  the  desired  parameters.  When 
executed,  the  user  is  prompted  to  input  the  number  of  simulation  runs  desired  and  to 
select  the  "actual  winds"  to  be  used.  The  program  then  randomly  selects  the  parameters 
identified  in  Table  4. 


Parameter 

Method 

Possible  Results 

Planned  Wind  File 

\ 

Uniform  Distribution 

any  valid  wind  file  from 
the  same  day  as  the 
selected  "actual  winds" 

SAON 

selects  which  sensor  model  is  used 

Uniform  Distribution 

YES/NO 

CARP  Offset 

distance  from  desired  release  point 

Uniform  Distribution 

0,  1000,  2000,  3000  feet 
each  axis 

PSIDOT 

parachute  turn  rate 

Normal  Distribution 

Mean:  1.89  degrees 

STDEV:  1  degree 

RELEASEALT 

altitude  at  parachute  release 

Based  on  maximum 

Height  of  actual  wind  file 

9500  feet,  20000  feet 

CompassBias 

Uniform  Distribution 

+  2  degrees 

SASEED 

seeds  used  in  GPS  error  models 

Uniform  Distribution 

0  to  9999 

Table  4.  Simulation  Parameter 


Another  script  file,  titled  CARLO  is  then  executed.  This  execution  file  imports 
the  parameters  in  the  initialization  file,  executes  the  simulation  CARP  to  calculate  the 
predicted  trajectory  and  the  simulation  C9POINTMASS  to  perform  the  simulation.  The 
results  are  then  stored  in  individual  data  files.  Summary  statistics  are  calculated  by  a 
script  titled  MONTE  called  from  within  the  CARLO  routine  and  stored  in  a  text  file. 
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C.  RESULTS 

Appendix  A  contains  the  initialized  states  and  results  for  each  of  the  600 
simulations  conducted.  These  data  include  the  achieved  accuracy  improvement  of  the 
controlled  system  over  that  of  an  uncontrolled  system.  The  simulations  produced 
excellent  results  with  an  accuracy  of  210  feet  (64  meters),  Circular  Error  Probable  (CEP). 
The  total  average  horizontal  error  was  309  feet  (94  meters)  with  an  average  of  15  control 
inputs  being  required.  The  maximum  number  of  control  inputs  for  all  600  simulations 
was  33.  Appendix  B  presents  sample  plots  from  several  simulations. 

Figures  52  through  54  present  three  dimensional  plots  of  several  of  the  simulation 
results.  These  figures  illustrate  the  initial  release  points  and  flight  paths  to  impact  with 
the  center  of  the  coordinate  system  (0,0)  being  the  planned  impact  point.  With  the 
exception  of  a  few  trials,  the  system  guided  to  the  desired  target. 


Figure  52.  Simulation  Results,  10000  foot  Release,  100  trials 
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Figure  53.  Simulation  Results,  20000  foot  Release,  100  trials 


Figure  54.  Simulation  Results,  50:  10000  foot  Release,  50:  20000  foot  Release 
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To  assess  each  of  the  randomized  variables,  the  statistics  were  summarized  by  the 
time  from  the  predicted  winds  to  the  actual  winds,  the  GPS  errors  introduced  (state  of 
selective  availability  (ON  =  included)),  and  the  distance  from  the  predicted  CARP  to 
simulated  release  point.  Table  5  presents  these  results. 


Release  Alt 

Age  of 

SA 

Number 

Mean  Horizontal  Error  (feet) 

Control  Inputs  | 

(feet) 

Wind  (hrs) 

On/Off 

Of  Rims 

Uncontrolled 

Controlled 

Average 

9500 

1 

0 

75 

1935 

159 

14 

3 

22 

9500 

1 

1 

91 

213 

14 

3 

21 

9500 

■  2 

0 

52 

2169 

157 

13 

3 

21 

9500 

2 

1 

50 

2232 

192 

13 

3 

18 

9500 

3 

0 

28 

2355 

56 

15 

3 

20 

9500 

3 

1 

28 

3340 

136 

16 

3 

22 

9500 

5 

0 

38 

3188 

401 

12 

5 

20 

9500 

5 

1 

51 

3412 

448 

11 

5 

20 

9500 

6 

0 

14 

3351 

427 

10 

3 

17 

9500 

6 

1 

22 

3816 

573 

10 

4 

15 

9500 

7 

o 

29 

4044 

818 

9 

3  • 

13 

9500 

1 

1 

22 

4326 

886 

10 

4 

16 

20000 

3 

0 

16 

3501 

82 

24 

4 

30 

20000 

3 

1 

15 

3372 

173 

22 

5 

31 

20000 

5 

0 

11 

3179 

312 

. TT 

5 

28 

20000 

5 

"1 

ib 

2923 

352 

21 

4 

25 

20000 

7 

0 

15 

5886 

138 

22 

6 

33 

20000 

7 

1 

13 

5999 

208 

23 

6 

33 

20000 

9 

0 

12 

7492 

828 

13 

5 

19 

20000 

9 

1 

8 

8081 

701 

18 

5 

24 

Table  5.  Simulation  Results 


The  predominant  factor  influencing  accuracy  of  the  control  system  was  the  time 
from  the  predicted  winds  used  to  establish  the  planned  trajectory  to  that  of  the  time  of  the 
simulated  airdrop.  Figure  55  plots  the  resultant  accuracy  with  the  “age”  of  wind  and 
shows  that  the  impact  accuracy  is  greatly  reduced  when  the  actual  winds  are  within 
approximately  three  hours  of  the  predicted  winds. 


83 


Figure  55.  Influence  of  Age  of  Wind  on  Accuracy 
Figure  56  presents  an  example  of  the  magnitude  of  wind  error  for  the  three  and 


five  hour  old  winds.  For  the  first  set  of  winds  (three  hours),  the  actual  winds  closely 
match  those  of  the  predicted  winds  (within  15  feet  per  second).  The  second  example 
(five  hours)  shows  differences  in  wind  speed  up  to  30-40  feet  per  second  with  changes  in 
significant  changes  in  direction  below  2,000  feet  in  altitude.  The  system  simply  does  not 


have  enough  drive  to  overcome  differences  of  this  magnitude  between  the  planned  winds 


and  the  actual  winds. 


Figure  56.  Wind  Comparison 
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The  effects  of  GPS  navigation  errors  are  not  as  apparent  as  that  of  age  of  wind. 
Figure  57  presents  the  accuracy  with  the  status  of  GPS  errors  (Selective  Availability). 
This  plot  shows  no  direct  correlation  between  the  state  of  GPS  errors  and  the  resultant 
impact  accuracy. 


Figure  57.  Influence  of  Selective  Availability  Errors  on  System  Accuracy 

Evaluating  the  trials  that  had  impact  errors  less  than  100  feet,  a  slight  correlation 
is  observed  between  Selective  Availability  status  and  impact  errors.  For  SA  Off  trials, 
the  resultant  impact  errors  had  a  smaller  variance  than  that  for  SA  On  (Figure  58).  This 
shows  that  the  GPS  navigation  errors  do  contribute  to  the  overall  accuracy  of  the  AGAS 
system  but  do  not  have  the  same  degree  of  impact  as  the  age  of  wind  data. 

The  effect  of  the  magnitude  of  the  offset  of  the  actual  release  position  from  that  of 
the  planned  release  point  also  does  not  show  a  direct  correlation  with  the  magnitude  of 
the  impact  accuracy  (Figures  58  and  59). 
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Figure  58.  Influence  of  Selective  Availability  Errors  on  System  Accuracy 


Figure  59.  Influence  of  CARP  Errors  on  System  Accuracy 


xn.  CONCLUSIONS  AND  RECOMMENDATIONS 

A.  CONCLUSIONS 

This  study  demonstrated  that  the  AGAS  concept  is  capable  of  providing  a  cost 
effective  option  for  precision  airdrop.  Simulations  demonstrate  that  accuracies  of  210 
feet  (64  meters).  Circular  Error  Probable  (CEP),  can  be  achieved.  The  development  is 
warranted  in  going  into  the  next  stage  of  development.  Additional  efforts  are  needed  to 
optimize  the  control  system  in  an  effort  to  further  reduce  the  amount  of  stored  gas 
required  on  the  operational  system. 

The  flight  test  program  provided  adequate  flight  dynamic  data  for  the  AGAS 
system.  The  instrumentation  system  developed  for  this  program  proved  invaluable.  The 
wind  measurement  techniques  employed  by  use  of  a  "calibration"  parachute  system 
provided  the  best  possible  wind  measurement  system  for  this  application.  Use  of  the 
measured  GPS  ground  track  velocities  as  the  wind  estimate  was  validated  to  provide 
sufficient  wind  measurements. 

The  MMLE3  parameter  estimation  algorithms  as  implemented  in  MATLAB386 
are  not  sufficient  for  predicting  the  stability  derivatives  for  this  parachute.  A  numerical 
sensitivity  to  the  magnitude  of  the  particular  parameters  is  suspected  to  be  the  cause.  It  is 
unlikely  that  these  techniques,  had  they  provided  a  solution,  would  have  produced  a 
reliable  solution  due  to  the  non-linearity  of  the  aerodynamics  of  a  parachute  system. 

An  efficient  Monte-Carlo  type  simulation  was  developed  using  a  point  mass 
model  for  parachute  dynamics,  sensor  models  for  GPS  and  heading  information,  and  a 
Bang-Bang  type  control  system. 
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The  point  mass  concept  for  system  modeling  is  adequate  for  this  evaluation.  The 
position  and  velocity  results  from  the  model  matched  well  with  the  measured  flight  test 
results  under  the  same  wind  conditions.  However,  to  optimize  the  control  system,  a  full 
six  degree  of  freedom  model  is  likely  required  to  capture  the  proper  heading  response.  A 
complete  set  of  general  6-DOF  equations  of  motion  were  developed  and  presented  herein. 
However,  additional  work  is  needed  to  either  estimate  the  stability  derivatives  for  the 
parachute  system  or  develop  a  physics-based  aerodynamic  model  such  as  that  proposed 
by  Colin  Tory.26 

Six  hundred  simulations  were  conducted  with  randomly  selected  initialization 
parameters.  These  results  demonstrated  that  the  Affordable  Guided  Airdrop  System,  as 
described  herein,  shows  strong  potential  of  providing  a  low-cost  alternative  for  precision 
airdrop.  Three  critical  factors  will  drive  the  final  design  of  the  AGAS.  First,  the 
accuracy  of  the  estimated  winds  when  determining  the  planned  trajectory  is  the  dominant 
factor  in  the  accuracy  of  the  AGAS  concept.  Winds  of  up  to  6  hours  old  (as  compared  to 
the  'actual'  winds  used  in  the  simulation)  resulted  in  large  horizontal  errors  from  the 
desired  impact  points.  Secondly,  the  rotation  rate  of  the  parachute  system  is  important. 
Rotation  rates  with  a  mean  of  1 .89  degrees  per  second  and  a  standard  deviation  of  1 
degree  per  second  allowed  effective  control.  If  the  rotation  rates  of  the  production 
system  are  increased  from  that,  sufficient  control  may  not  be  possible.  Finally,  the 
number  of  control  inputs  required  to  achieve  the  desired  accuracy  is  marginal  under  the 
current  control  concept.  However,  no  attempts  were  made  to  optimize  the  control 
algorithms  for  minimum  fuel  usage. 
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B.  RECOMMENDATIONS 


This  study  demonstrated  the  feasibility  of  the  AGAS  concept  to  provide  a  low 
cost  alternative  for  precision  airdrop.  However,  the  success  of  the  final  design  rides  on 
the  three  critical  factors  presented  above.  Therefore,  the  following  recommendations  are 
offered  for  follow-on  work: 

1.  Complete  the  incorporation  of  6-DOF  equations  of  motion.  Nonlinear 
parameter  estimation  techniques  would  need  to  be  investigated  or  a  physics-based 
approach  would  be  needed. 

2.  Fully  characterize  the  performance  of  the  AGAS  concept  using  the  G-12 
and/or  G-l  1  parachute  systems.  The  remote  controlled  activation  technique  used  on  the 
C-9  test  program  should  be  applied  to  the  G-12  system  as  soon  as  G-12  actuators  are 
available. 

3.  Investigate  optimizing  the  control  algorithms  for  minimum  fuel  usage.  The 
current  methodology  provides  minimum  horizontal  errors  without  regard  to  fuel 
consumption.  The  fuel  consumption  for  these  algorithms  is  marginal.  A  technique  for 
evaluating  multiple  control  algorithms  is  the  application  of  a  ground  based  guidance 
computer.  The  navigation  data  from  the  system  could  be  downlinked  via  radio  modem  to 
the  guidance  computer,  the  control  algorithms  could  determine  the  desired  activation  of 
the  actuator,  and  these  control  commands  could  be  uplinked  to  the  test  item.  Figure  60 
illustrates  this  concept. 
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Figure  60.  Autonomous  Guidance  Test  Concept 
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APPENDIX  A.  SIMULATION  RESULTS 
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APPENDIX  B.  SAMPLE  SIMULATION  RESULTS 


This  appendix  contains  output  for  numerous  simulation  trials  to  illustrate  a  sample 
of  the  results  obtained. 


Run35 

X  Error:  -2.8414  Y  Error:  -11.3694  Hor  Error  11.7191 
No  Control  X  Error:  -466.9893  Y  Error  -122.2171  Hor  Error  482.7174 
A  Inputs:  4  B  Inputs:  3  C  Inputs:  6  D  Inputs:  2 
Total  Control  Inputs:  15 


Run202 

X  Error:  688.6045  Y  Error:  -2295.7423  Hor  Error  2396.7913 
No  Control  X  Error:  5963.0562  Y  Error:  -541 1.0036  Hor  Error:  8052.1425 
A  Inputs:  1  B  Inputs:  1  C  Inputs:  0  D  Inputs:  1 
Total  Control  Inputs:  3 


Run217 

X  Error: -209.5789  Y  Error  -1305.1452  Hor  Error:  1321.8651 
No  Control  X  Error:  7971.5228  Y  Error:  -7412. 1 185  Hor  Error:  10885.0667 
A  Inputs:  6  B  Inputs:  5  C  Inputs:  5  D  Inputs:  6 
Total  Control  Inputs:  22 


Run217  (continued) 


Control  Activation:  Run  217 


Run223 

X Error: -115.2622  YError:  -37.2997  HorError:  121.1473 
No  Control X Error: -684.7159  YError:  1083.0916  HorError:  1281.3755 
A  Inputs:  6  B  Inputs:  9  C  Inputs:  7  D  Inputs:  8 
Total  Control  Inputs:  30 


iiiiliii 
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Run223  (continued) 
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APPENDIX  C.  SIMULINK®  REALIZATION 


The  following  presents  the  Simulink®  realization  for  this  effort.  It  is  presented  in 
the  hierarchy  of  the  implementation  of  the  model. 


MUST  LOAD  pmdata.mat  BEFORE  RUNNING 


Top-level  Model 


Dynamic  Model 
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CompassModel 


Overall  Sensor  Model 


GPS  Error  Models 


NOTE:  Trigger  Controls  Timing  to  Permit  Control  Changes 


Controller 
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Controller  ‘Which  Control’  Block:  Determines  Which  Actuator  to  Affect 


Controller  ‘On_OfF  Block:  Checks  Tolerances  and  Activates  Control 
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APPENDIX  D.  AGAS  INSTRUMENTATION  DESCRIPTION 


Overview 

The  AGAS  instrumentation  package  was  required  to 
measure  and  record  the  state  of  the  AGAS  payload  to 
include:  position,  velocity,  acceleration,  and  three- 
axis  attitude  and  attitude  rates.  In  addition,  die  state 
(extended  or  retracted)  of  each  of  the  four  actuators 
was  recorded.  The  position  and  velocity  of  the  AGAS 
payload  was  provided  by  a  differential  carrier-phase 
GPS  system  known  as  the  Improved  Vehicle 
Tracking  System  (IVTS).  Acceleration  data  were 
derived  from  a  triad  of  accelerometers,  and  the 
attitude  of  the  package  was  measured  by  an  Attitude 
Heading  Reference  System  (AHRS).  Pressure 
transducers  attached  to  the  PMAs  provided  the  state 
of  the  actuators.  Since  time  correlation  of  the  data 
from  multiple  sources  was  critical,  a  timing  source 
was  included  in  the  package.  All  measured  data  were 
recorded  on  removable  flash  memory  devices. 

The  AGAS  instrumentation  package  contained  two 
PC/104  computer  systems.  The  sensor  computer 
system  managed  the  data  acquisition  and  recording 
for  the  AHRS,  acceleration  and  actuator  states 
collected  at  a  15-Hertz  (Hz)  rate.  The  Global 
Positioning  System  (GPS)  computer  system 
controlled  the  GPS  receiver  and  acquired  and 
recorded  the  receiver  measurements.  The  GPS 
trajectory  solution  was  recorded  at  a  5-Hz  rate. 
Details  of  the  instrumentation  system  design  are 
presented  in  the  appendix.  The  system  design  is 
illustrated  in  Figure  1. 


Figure  1.  Instrumentation  Overview 


System  Architecture 

The  data  acquisition  system  was  based  on  an  industry 
standard  embedded  computer  product  family  called 
PC/104.  This  product  family  consists  of  the  central 
processing  units,  power  supplies,  and  video  as  well  as 


a  large  assortment  of  analog  and  digital  data 
acquisition  and  control  modules  useful  in  configuring 
data  acquisition  systems.  The  CPU  modules  are 
based  on  Intel®  X86  processors  and  use  the  Industry 
Standard  Architecture  (ISA)  bus  associated  with  the 
IBM®  PC/AT  series  of  personal  computers.  Typical 
enhancements  to  the  standard  computer  include  solid 
state  mass  storage,  Basic  Input/Output  Systems 
(BIOS)  extensions,  and  reduced  power  requirements. 
These  modules  are  highly  integrated  and  have  a  small 
footprint,  measuring  3.6  X  3.8  inches.  The  stacking 


Figure  2.  PC-104  Modules 


feature  of  the  PC/104  family  makes  it  easy  to 
configure  a  system  to  meet  unique  requirements.  A 
typical  PC/104  system  is  shown  in  Figure  2. 

AGAS  Instrumentation  Computer  System 

The  sensor  computer  system  contained  a  central 
processing  unit  (CPU)  module,  an  analog-to-digital 
converter  (ADC)  module,  a  power  supply  module,  a 
timing  module  and  a  type  II  PCMCIA  module.  The 
CPU  module,  is  a  highly  integrated  module 
containing  an  Intel  486  DX4  processor  operating  at 
100  Mhz,  an  IDC  hard  disk  controller,  a  floppy  disk 
controller,  two  serial  ports,  a  parallel  port  and  4 
megabytes  of  random  access  memory.  The  module 
also  contained  a  1  megabyte  flash  memory  which, 
through  BIOS  extensions,  looked  to  the  processor 
like  a  bootable  IDC  hard  disk.  This  flash  memory 
was  programmed  with  the  operating  system  and  the 
application  software  that  controlled  the  modules  on 
the  stack. 

The  ADC  module  is  an  8  channel  12-bit  converter 
configured  for  5-volt  bipolar  input  An  internal 
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adjustable  gain  amplifier  provides  programmable 
gains  of  1,  10,  100,  or  1000.  The  analog  to  digital 
conversion  rate  is  a  function  of  the  mode  used  to 
initiate  the  conversion,  and  varies  from  2,000 
conversions  per  second  in  the  software  conversion 
mode  to  a  maximum  of,  100,000  conversions  per 
second  if  the  direct  memory  access  (DMA) 
conversion  mode  is  used.  An  interrupt  conversion 
mode  is  also  available,  which  will  provide  up  to 
20,000  conversions  per  second.  In  this  mode,  the 
interrupt  signal  may  be  provided  by  an  external 
signal  or  it  may  be  provided  from  the  output  of  a 
programmable  32-bit  counter  contained  on  the 
module.  This  module  also  contains  two  12-bit 
digital-to-analog  converters  as  well  as  four  digital 
inputs  and  four  digital  outputs. 

The  timing  module  is  a  digital  clock  that  provides 
accurate  time  to  application  software  through  the 
PC/104  bus.  The  clock  may  also  be  initialized  from 
the  application  software  through  the  bus.  Time-of- 
year  is  available  in  a  binary  coded  decimal  (BCD) 
format  with  a  resolution  of  hundreds  of  nanoseconds. 
Clock  accuracy  is  maintained  by  phase-locking  the 
internal  clock  to  an  external  reference.  The  reference 
may  be  any  of  the  standard  timing  signals  usually 
found  on  United  States  test  ranges,  such  as  the  Inter- 
Range  Instrumentation  Group  (IRIG)  time  code 
signals  or  NASA  time  code  signals  in  addition  the 
internal  clock  may  be  synchronized  to  an  external 
one  pulse  per  second  (1-PPS).  Synchronization 
accuracy  to  the  IRIG  or  NASA  timing  signals  is 
specified  to  be  between  20  microseconds  and  5 
microseconds  depending  upon  which  code  is  being 
used  for  the  synchronization  input.  A 
synchronization  accuracy  of  one  microsecond  is 
specified  for  the  1-PPS  input. 

The  sensor  computer  system  power  supply  module  is 
a  DC  to  DC  converter  that  provides  up  to  50  watts  of 
power  in  the  form  of  +5  volts  and  +-  12  volts.  The 
power  supply  module  requires  a  DC  input  voltage 
between  6  volts  and  40  volts.  Each  of  the  power 
signals  is  connected  to  the  PC/104  bus  to  facilitate 
powering  the  modules  in  the  stack. 

The  PCMCIA  module  located  in  the  sensor  computer 
module  stack  is  a  standard  type  II  PCMCIA  interface. 
With  the  appropriate  driver  software,  flash  memory 
modules  inserted  into  the  interface  appear  to  the 
operating  system  as  an  IDC  hard  drive. 

Three  modules  comprised  the  GPS  system  computer: 
a  CPU  module,  a  power  supply  module,  and  a 
PCMCIA  module.  The  hardware  on  these  three 
modules  is  identical  to  the  hardware  in  the  sensor 


system  computer  modules.  The  CPU  was 
programmed  with  the  software  required  to  interface 
with  the  GPS  receiver  sensor.  A  photograph  of  the 
packaged  instrumentation  system  is  included  in 
Figure  3. 


Figure  3.  AGAS  Instrumentation  Package 


AGAS  Sensors 

As  mentioned,  the  AGAS  instrumentation  package 
contained  four  types  of  sensors:  a  GPS  receiver,  an 
AHRS  unit,  an  accelerometer  triad,  and  muscle  state 
sensors.  The  GPS  sensor  consists  of  an  LI  C/A  code 
receiver  and  antenna.  The  receiver  contains  12 
independent  channels,  which  track  the  GPS  satellite 
signal  in  parallel.  The  receiver  is  capable  of  making 
measurements  at  a  20  Hz  rate  and  providing  position 
and  velocity  output  up  to  a  10  Hz  rate  (a  5  Hz  rate 
was  used  for  the  AGAS  implementation).  Two  RS- 
232C  I/O  ports  are  available  on  the  receiver  to 
provide  for  control  inputs  and  data  outputs.  Either  of 
the  I/O  ports  may  be  used  to  input  commands  and 
output  data,  or  both  may  be  used  at  the  same  time.  An 
auxiliary  I/O  port  provides  a  1PPS. 

The  AHRS  used  in  the  AGAS  instrumentation  system 
is  a  solid  state  device  that  measures  three-axis 
attitude  rate  and  integrates  it  to  form  attitude  and 
heading.  The  attitude  and  heading  is  compared  with 
two  reference  vertical  pendulums  and  a  tri-axial  flux- 
gate  magnetometer.  The  resulting  error  is  filtered 
and  used  to  adjust  the  output  of  the  system  causing  it 
to  converge  to  the  attitude  of  the  vertical  pendulums 
and  to  magnetic  heading.  The  attitude  and  attitude 
rates  are  available  on  an  RS-232C  interface  at  a  rate 
of  just  less  than  15  Hz.  The  specified  accuracy  of  the 
AHRS  data  is  shown  in  Table  1. 
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Static 

Dynamic 

(percent) 

Rate  Accuracy  (deg/sec) 

+  0.2 

+  2 

Attitude  Accuracy  (deg) 

+  0.5 

+  2 

Heading  Accuracy  (deg) 

+  1.0 

+  2 

Table  1.  AHRS  Accuracies 


The  accelerometers  in  the  AGAS  instrumentation 
system  are  miniature  devices  with  a  full-scale  output 
of  +  lOg.  The  acceleration  experienced  by  these 
devices  is  converted  to  an  electrical  signal  via  strain 
gauges  wired  as  a  Wheatstone  bridge  and  attached  to 
the  acceleration-sensing  element.  The  frequency 
response  of  these  devices  is  nominally  0  to  300  Hz, 
with  a  minimum  response  of  140  Hz.  The  sensitivity 
of  each  element  is  approximately  8  milli-volts/g. 
Three  accelerometers  were  mounted  in  a  steel 
mounting  block  with  their  input  axis  orthogonal.  One 
axis  was  defined  to  be  the  X-axis  with  the  +  input  of 
the  accelerometer  in  the  direction  of  +X.  This  axis 
was  designated  as  the  axis  pointing  forward.  The 
other  two  devices  were  mounted  such  that  the  +Y 
axis  was  pointing  to  the  right  and  the  +Z  axis  was 
pointing  down.  The  mounting  block  was  attached  to 
the  mounting  structure  with  the  +X  axis  of  the 
accelerometer  triad  parallel  to  the  AHRS  +X  axis. 

The  last  group  of  sensors  in  the  AGAS  sensor  system 
is  the  riser  muscle  state  sensors.  These  were  simply 
four  pressure  transducers,  one  for  each  muscle.  The 
four  muscles  in  the  control  system  have  two  states, 
extended  and  contracted.  When  the  muscle  is 
contracted,  the  pressure  in  the  lines  supplying  air  to 
the  muscles  is  about  50  psi.  When  the  muscle  is 
extended,  the  pressure  drops  to  nearly  zero. 
Therefore,  the  pressure  in  the  lines  represents  the 
state  of  the  muscles  (zero  indicating  the  control  is 
activated). 

AGAS  Instrumentation  Software 

The  software  residing  in  the  processors  shown  in  the 
diagram  controlled  the  hardware  in  the 
instrumentation  package.  Both  processors  ran  the 
MSDOS®  operating  system,  and  the  application 
software  was  written  in  the  C  language. 

The  application  software  in  the  GPS  processor  was 
the  same  as  that  used  in  the  Yuma  Proving  Ground’s 
Improved  Vehicle  Tracking  System,  and  in 
conjunction  with  the  receiver  used,  will  produce  sub¬ 
meter  positioning  when  the  data  is  post-processed  to 
remove  errors.  This  software  initialized  the  GPS 
receiver  at  power  on  and  began  to  acquire  the  GPS 
satellite  signal.  While  the  AGAS  package  was  in  the 
aircraft,  the  GPS  receiver  tracked  the  signal  via  a 


retransmission  system  that  received  the  signal 
external  to  the  aircraft  and  retransmitted  it  inside  the 
aircraft.  Once  the  GPS  receiver  was  tracking 
satellites,  the  processor  collected  the  navigation  data 
and  recorded  it  on  the  recorder. 

The  sensor  instrumentation  processor  also  began 
initializing  at  power  up.  Early  in  this  process,  it 
started  requesting  a  time  message  from  the  GPS 
receiver.  Once  this  message  was  received  with  status 
indicating  the  GPS  receiver  had  solved  for  GPS  time, 
the  sensor  instrumentation  processor  initialized  the 
timing  module  to  the  current  UTC  time.  The 
accuracy  of  the  time  in  the  timing  module  was 
subsequently  maintained  by  phase  locking  its  clock  to 
a  1 PPS  signal  from  the  GPS  receiver.  Once  accurate 
time  had  been  established  in  the  timing  module,  the 
sensor  processor  enabled  the  AHRS  output  and  began 
acquiring  data. 

The  attitude  and  attitude  rate  data  from  the  AHRS 
was  transferred  to  the  processor  through  an  RS-232C 
interface.  Each  transfer  contained  a  string  of  ASCII 
characters  that  included  a  unique  header,  the  current 
AHRS  attitude  and  attitude  rate.  The  transfer  was 
terminated  with  the  ASCII  character  representing  a 
carriage  return.  A  period  of  inactivity  on  the 
interface  separated  the  transfers.  At  the  point  the 
software  detected  the  unique  character  designating 
the  start  of  a  transfer,  the  current  time  was  sampled 
and  stored.  This  time  became  the  time  of  the  attitude 
sample. 

As  soon  as  the  attitude  time  sample  was  stored,  the 
software  acquired  the  accelerometer  and  pressure 
data  from  the  ADC.  This  activity  started  with 
another  sample  of  current  time,  which  became  the 
sample  times  of  the  data  from  the  A/D  converter. 
After  the  accelerometer  and  pressure  data  were 
sampled,  the  software  converted  the  data  to  ASCII 
and  recorded  it  in  a  file  designated  for  A/D  converter 
data. 

At  this  point  the  software  continued  with  the 
acquisition  of  the  attitude  record.  Since  the  serial 
port  driver  •  operates  under  interrupt  control,  the 
attitude  data  that  arrived  at  the  port  during  the  A/D 
sample  interval  was  placed  in  a  buffer.  These  data 
were  retrieved,  formatted,  and  recorded  in  the  file 
designated  for  attitude  data.  This  finished  one  loop 
of  the  data  acquisition  process  and  the  software 
began  to  monitor  the  serial  port  for  the  beginning  of  a 
new  attitude  sample.  The  message  rate  from  the 
AHRS  was  approximately  15  Hz,  so  the  sample  rate 
for  all  data  was  also  at  approximately  15  Hz. 
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APPENDIX  E.  MATLAB®  SCRIPT  FILES 

The  following  presents  the  MATLAB®  script  files  utilized  in  this  effort.  Note 
that  only  the  critical  software  code  is  presented. 

Initialization  File 

%This  program  initializes  the  workspace  for  running  the 
simulation. 

load  pmdata 

runs=input( ’Desired  Number  of  Runs:  '); 

disp ( 'To  accept  current  parameters  (one  by  one)  just  hit 
return ' ) ; disp ( '  ’ ) ; disp ( '  ' ) ; 

disp ('Current  Windfile:  ' ) ;disp (windfile) ; 

newwindfile=input (' Select  a  Wind  File  ([day  drop#]  eg.  [310  1]): 
newwindfile 

if  isempty (newwindfile) ~=1 
windf ile=newwindf ile ; 
end 

for  i=l:runs 

if  windfile (1,1) ==310; 

RAND=ceil ( rand ( 1 ) *3.999) ; 

filechoice=[ 110607  110609  110612  110613]; 

end 

if  windfile  (1, 1)  =311; 

RAND=ceil ( rand ( 1 ) *2.999)  ; 
f ilechoice= [110707  110711  110714]; 

end 

if  windfile (1, 1) ==0617; 

RAND=ceil (rand(l) *3.999)  ; 
filechoice= [ 061710  061712  061714  061716]; 
end  _ _ _ 
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rawinfile=filechoice(RAND) ; 
saon=floor  (rand (1)  *1.999)  ,- 
%sastart=ceil (rand (1) *6. 999) ,- 
RAND=ceil (rand  (1) *3 .999)  ; 

offset*  [0  0,-1000  1000; -2000  2000,-3000  -3000],- 
xof fset=off set  (RAND,  1)  ,-yof f set=of f set  (RAND,  2)  ; 
psidot=randn(l)+l.85; 

BiasDir=floor  (rand  (1)  *1.99)  ,- 
if  BiasDir==l 
CompassSign=l; 
else 

Compass  S ign= - 1 ,- 

end 

CompassBias* (rand(l) *1.999) *CompassSign; 
if  windf ile (1, 1) ==617; 

releasealt=20000; 

else 

releasealt =95  00,- 

end 

if  exist (' initf ile ') ==0 

initfile= [windf ile  rawinfile  saon  releasealt  Vi  Hi  dzalt 
xoffset  yoffset  psidot  CompassBias]  ,- 
save  init  initfile 
else 

ile  (i ,;)  =  [windf  ile  rawinfile  saon  releasealt  Vi  Hi 

- dzalt  xoffset  yoffset  psidot  CompassBias]  ,- 

save  init  initfile  -append 

end 

clear  new*  acchange 
end 

disp ( ' init  done ' ) 


Equations  of  Motion  -  dof3.m 


function  accel=dof3(z); 

V=z(l:3); 

Vt=norm(V); 

%checks  to  see  if  density  is  contained  in  the  input  vector. 

%If  not,  sets  to  standard  sea  level  density 
if  length(z)==5; 
rho=z(4); 
time=z(5); 
else 

rho=0.002377; 

time=z(5); 

end 

%Atmospherics 

qbar=.5*rho*(VtA2); 

g=32.17; 

%System  aerodynamics/characteristics 

opentime=3;  %seconds 

Cd=0.75;  %obtained  from  flight  test 

Do=28;  %reference  diameter  in  feet 

Dp=.67*Do;  %profile  diameter 

So=pi/4*DoA2;  %reference  area  (flat  circular  parachute) 

W=346; 

m=W/g; 

%Calculate  drag  area.  Assumes  a  linear  opening  of  the  parachute 
if  time<=opentime 
CdSt=[0  0;opentime  Cd*So]; 
CdSo=interpl(CdSt(:,l),CdSt(:,2),time); 
else 

CdSo=Cd*So; 

end 
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%Calculate  mass  terms 

alpha  1 1 =0.25*rho*4/3*pi*(Dp/2)  A3; 

alpha22=alphall; 

alpha33=2*alphall; 

Ml=[(m+alphall)  0  0; 

0  (m+alpha22)  0;  • 

0  0  (ra+alpha33) ] ; 

%accel= [-qbar*CdSo/m*V ( 1) /Vt ; -qbar*CdSo/m*V (2) /Vt ; - 
qbar*CdSo/m*V(3)/Vt+g] ; 

acce 1 = inv (Ml ) * ( [ - qbar* CdSo*V ( 1 ) /Vt ; - qbar* CdSo*V ( 2 )  /Vt 
-qbar*CdSo*V(3)/Vt]+[0;0;W] ) ; 


Monte-Carlo  Simulation  carlo.m 

%This  file  executes  the  desired  simulations 
clear 

load  pmdata;load  init;[row  col]=size(initfile); 

startrun=input('Key-in  Starting  Run  Number  ’); 
for  runno=l:row 
windfile=initfile(runno,l:2); 

eval([’Wind=Wind'  num2str(windfile(l))  'd'  num2str(windfile(2)) 

rawinfile=initfile(runno,3); 

eval([’rawin=rawin'  num2str(rawinfile) 

saon=initfile(runno,4); 

if  saon==l 

saseed=[floor(rand(l,3)*1000)zeros(l,6)]; 

else 

saseed=floor(rand(l,9)*1500); 

end 


releasealt=initfile(runno,5); 

Vi=initfile(runno,6); 

Hi=initfile(runno,7); 

dzalt=initfile(runno,8); 

xoffset=initfile(runno,9); 

yoffiset=initfile(runno,10); 

psidot=initfile(runno,ll); 

CompassBias=initfile(runno,12); 

sim('carp') 

disp([num2str(runno) '  CARP  Sim(s)  Done']); 
CARP=-carptraj(length(carptraj),l  :2); 
fliptraj 

sim('c9pointmass') 

%noisel=wavread('c:\windows\media\ofRce97\gunshot.wav'); 

%sound(noisel) 

disp([num2str(runno) '  Control  Sim(s)  Done']) 
monte; 

clear  carptraj; 
horimpacterr 

results(runno,:)=[runno+startrun-l  windfile  rawinfile  saon  ... 
xoffeet  yoffset  psidot  releasealt  Vi  Hi  dzalt  CARP  ... 
hornocontimpacterr  horimpacterr  totcontrols]; 
end 

eval(['save  runfile'  num2str(startrun) '_'  num2str(startrun+runno-l)... 
'.asc  results  -ascii']) 
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